這個案例由於原始數據過大，沒辦法進行原始數據預處理的朋友可以在limma包下載頁面下載已經預處理的數據，進行從15.1.6開始的統計分析。
檢測來源：不同的69約魯巴人的B淋巴細胞的細胞組，這個文庫是國際人體基因組計劃的一部分。
每個個體的RNA至少在Illumina Genome Analyser 2的兩個pannel上測量。在這個案例裡我們要尋找男性和女性的差異表達基因。這個案例研究需要limma3.9.19或以上版本。
15.1.2 數據獲取已預處理的數據可以在http://bioinf.wehi.edu.au/limma/ 下載。未處理的數據FASTQ 文件可以在下面的地址下載。
15.1.3 約魯巴個體和FASTQ文件：Pritchard Lab’s eQTL序列信息地址：http://eqtl.uchicago.edu/RNA_Seq_data/unmapped_reads/
樣本信息地址：http://eqtl.uchicago.edu/RNA_Seqdata/list_lanes_pickrell_2010_plosgenetics
性別信息（International HapMap Project NHGRI Repository ）：http://ccr.coriell.org/
我們只使用在Argonne國家實驗室的測序數據。Argonne中心的所有跑道都用46個bp讀數進行了排序。
> Targets <- read.delim("targets.txt", sep = "\t")> Targets center flow_cell lane individual concentration cellline gender1 argonne 090617_HGAC_S100091 5 NA18486 3.5 NA18486 Male2 argonne 090330_HGAC_S100065 3 NA18498_2 3.5 NA18498 Male3 argonne 090218_HGAC_S100047_PIPELINE_V2 7 NA18498 2.5 NA18498 Male#共有86個樣本，這裡只列出3個作為參考。15.1.4 將reads與基因組對應需要下載的文件：
基因組：GRCh37 (hg19 Genome Reference Consortium Human Build 37)
read 注釋文件：NCBI RefSeq Build 37.2. FASTA files。從UCSC Human Genome 的網址下載：http://hgdownload.soe.ucsc.edu/goldenPath/hg19/
> library(Rsubread)> FASTQFiles <- paste0(Targets$individual, "_", Targets$center, ".fastq.gz")> BAMFiles <- paste0(Targets$individual, "_", Targets$center, ".subread.bam")> buildindex(basename="hg19_subread",reference="GRCh37.fa")> align("hg19_subread", readfile1 = FASTQFiles, input_format = "gzFASTQ", output_file = BAMFiles)> gene <- featureCounts(BAMFiles, useMetaFeatures=TRUE, annot.inbuilt="hg19", allowMultiOverlap=TRUE)> propmap <- propmapped(BAMFiles)
##顯示對得上的reads數目，它的平均佔比是0.786。樣本用測序的cDNA的濃度進行排序。> summary(propmap[,"PropMapped"])Min. 1st Qu. Median Mean 3rd Qu. Max.0.597 0.778 0.798 0.786 0.810 0.828
> Targets$concentration <- factor(Targets$concentration)> o <- order(Targets$concentration)> barplot(as.matrix(t(cbind(propmap[o,"NumMapped"]*1e-6,+ (propmap[o,"NumTotal"]-propmap[o,"NumMapped"])*1e-6))),+ ylab="Number of Reads (millions)", names = Targets$individual[o],+ las=3, cex.names=0.3, cex.axis=0.7, density = c(80,0),+ legend = c("Number of Mapped Reads", "Number of Unmapped Reads"),+ args.legend = list(x = "topright",cex = 0.5, bty = "n"))加載基因水平計數
> Counts <- gene$counts> dim(Counts)[1] 25702 86> Counts[1:5,1:5] NA18486 NA18498_2 NA18498 NA18499 NA18501653635 12 8 19 17 5100422834 0 0 0 0 0645520 0 0 0 0 079501 0 0 0 0 0729737 6 6 2 6 4文庫的大小從60多萬到900多萬不等。樣品按用於測序的cDNA (pM)濃度排序。
> col.concent <- c("gold" , "orange", "red")> barplot(colSums(Counts[,o])*1e-6, names = Targets$individual[o], ylab="Library size (millions)",+ col = col.concent[Targets$concentration][o], las = 2, cex.names = 0.3)> legend("topleft", legend = c("1.5","2.5","3.5"), col = col.concent, pch = 15, cex = 0.7)為了在RNA-seq數據上運用limma，我們將比較女性和男性個體。
> table(Targets$gender[!duplicated(Targets$cellline)])Female Male41 3315.1.5 注釋注釋的Entrez Gene ID來自NCBI gene info file.
> library(limma)> columns <- c("GeneID","Symbol", "chromosome", "type_of_gene")> GeneInfo <- read.columns("130811-Homo_sapiens.gene_info",required=columns, stringsAsFactors = FALSE)> m <- match(gene$annotation$GeneID, GeneInfo$GeneID)> Ann <- cbind(gene$annotation[, c("GeneID", "Chr", "Length")],+ GeneInfo[m, c("Symbol", "type_of_gene")])> Ann$Chr <- unlist(lapply(strsplit(Ann$Chr, ";"), function(x) paste(unique(x), collapse = "|")))> Ann$Chr <- gsub("chr", "", Ann$Chr)> head(Ann) GeneID Chr Length Symbol type_of_gene28373 653635 1 1769 WASH7P pseudo37837 100422834 1 138 MIR1302-10 miscRNA27617 645520 1 1130 FAM138A miscRNA14957 79501 1 918 OR4F5 protein-coding29525 729737 1 3402 LOC729737 miscRNANA 100507658 1 793 <NA> <NA>15.1.6 DGEList 對象創建DGEList對象，把計數和相關的注釋結合起來：
> library(edgeR)> y <- DGEList(counts = Counts, genes = Ann)15.1.7 過濾##保留總計數大於50的基因：> A <- rowSums(y$counts)> isexpr <- A > 50##保留有注釋信息的基因> hasannot <- rowSums(is.na(y$genes)) == 0> y <- y[isexpr & hasannot, , keep.lib.size = FALSE]> dim(y)[1] 16918 8615.1.8 Scale標準化：> y <- calcNormFactors(y)##畫下MDS圖，圖像上更清晰地看到男性和女性的文庫區別> Gender <- substring(Targets$gender,1,1)> plotMDS(y, labels=Gender, top=50, col=ifelse(Gender=="M","blue","red"), gene.selection="common",+ prior.count = 5)從MDS圖中可以明顯看出低cDNA濃度和高cDNA濃度的分離。
> plotMDS(y, labels=Gender, col=col.concent[Targets$concentration], prior.count=5)> legend("topleft", legend = c("1.5","2.5","3.5"), col = col.concent, pch = 15)15.1.9 線性模型##設計設計矩陣> design <- model.matrix(~ gender + concentration, data = Targets)##用voom（）函數將read讀書轉換為log2-cpm，以及相應的權重，為後續線性模型服務。> v <- voom(y, design)##估計重複的細胞系之間的相關性，這些細胞系被多次測序> cor <- duplicateCorrelation(v, design, block = Targets$cellline)> cor$consensus[1] 0.446##細胞系內的相關性會略微改變voom的權重，所以我們再次運行voom:> v <- voom(y, design, plot = TRUE, block = Targets$cellline, correlation = cor$consensus)
同樣，我們更新voom權重的相關性
> cor <- duplicateCorrelation(v, design, block = Targets$cellline)> cor$consensus[1] 0.446可以看到，在第二次迭代中，相關性幾乎沒有改變。
現在尋找基因在男性和女性之間的差異表達。正對數變化意味著男性更高。
> fit <- lmFit(v, design, block = Targets$cellline, correlation = cor$consensus)> fit <- eBayes(fit)> summary(decideTests(fit))(Intercept) genderMale concentration2.5 concentration3.5Down 1464 27 1469 4105NotSig 3352 16869 13622 9704Up 12102 22 1827 3109> topTable(fit, coef = "genderMale", n = Inf, sort = "p", p = 0.05)[,-1]Chr Length Symbol type_of_gene logFC AveExpr t P.Value adj.P.Val B64595 Y 5242 TTTY15 miscRNA 6.529 -0.258 52.97 3.63e-68 6.14e-64 123.9836192 Y 897 RPS4Y1 protein-coding 9.691 2.224 44.08 1.91e-61 1.62e-57 117.2159086 Y 1390 EIF1AY protein-coding 7.962 0.659 42.44 4.52e-60 2.55e-56 113.6618284 Y 5581 KDM5D protein-coding 7.148 2.566 41.47 3.08e-59 1.30e-55 115.8028653 Y 4795 DDX3Y protein-coding 6.439 3.547 39.30 2.68e-57 9.08e-54 114.1557503 X 19271 XIST miscRNA -10.223 4.135 -37.53 1.20e-55 3.37e-52 107.001##結果過多，不一一列舉和預期的一樣，高表達的基因存在於性染色體上：
> chrom <- fit$genes$Chr> plotMD(fit, column=2, status=chrom, values=c("X","Y", "X|Y"),+ hl.col=c("red","blue", "green3"), main="Male vs Female",legend="bottomright")> abline(h=0,col="darkgrey")15.1.10 基因子集檢測獲取Y染色體男性特異性區域的基因列表，以及X染色體上據報導可以逃脫X抑制活性的基因列表。性染色體基因文件網址：http://bioinf.wehi.edu.au/software/GenderGenes/. 我們認為第一個列表中的基因在男性中應該是上調的，而第二個列表中的基因在女性中應該是上調的。
> load(url("http://bioinf.wehi.edu.au/software/GenderGenes/GenderGenes.RData"))> Ymale <- v$genes$GeneID %in% msYgenes> Xescape <- v$genes$GeneID %in% XiEgenesRoast基因組測試證實，在我們對男性和女性的比較中，男性特異性基因作為一個組顯著上升，而X基因作為一個組顯著下降：
> index <- list(Y = Ymale, X = Xescape)> mroast(v,index,design,2, block = Targets$cellline, correlation = cor$consensus, nrot=9999) NGenes PropDown PropUp Direction PValue FDR PValue.Mixed FDR.MixedY 13 0.000 1.0000 Up 1e-04 1e-04 1e-04 1e-04X 54 0.574 0.0926 Down 1e-04 1e-04 1e-04 1e-04
來自競爭性camera基因組測試的結果甚至更有說服力：
> camera(v,index,design,2)NGenes Direction PValue FDRY 13 Up 2.40e-262 4.80e-262X 54 Down 5.58e-26 5.58e-26我們可以在均值差圖上突出男性特異性基因和x逃逸基因：
> status <- rep("Other",nrow(fit))> status[Ymale] <- "Ymale"> status[Xescape] <- "Xescape"> plotMD(fit,column=2,status=status,values=c("Xescape","Ymale"),hl.col=c("red","blue"),+ main="Male vs Female",legend="bottomright")> abline(h=0,col="darkgrey")另一種查看基因的方式是查看Y和X特定基因在差異表達結果中的排名。這表明，除了少數例外，Y特異性基因的排位較高(陽性)，而X特異性基因的排位較低(陰性)：
> barcodeplot(fit$t[,2],Ymale,Xescape, labels=c("Up in male","Down in male"))> legend("top",legend=c("Ymale genes","Xescape genes"),lty=1,col=c("red","blue"))
15.2 Pasilla敲除後的不同剪接15.2.1這個案例研究重新分析了Brooks等人的數據。NOVA1和NOVA2蛋白是已知的調節哺乳動物剪接的蛋白。Brooks等以melanogaster果蠅為模型系統，利用RNA幹擾系統(RNA interference system, RNAi)敲除NOVA1和NOVA2的D. melanogaster or同源物表達，即Pasilla (ps)。實驗比較了經過處理和未經處理的S2-DRSC細胞系細胞。
在本例研究中，我們發現了在ps敲除後表達差異的外顯子和基因，我們發現敲除後的樣品與野生型相比，存在基因剪接差異。
為了完整起見，我們將描述從原始數據文件進行分析的所有步驟，包括映射和對齊。如果您想跳過映射並讀數步驟，請直接進入18.2.6節，即統計分析開始的地方。
15.2.2 GEO樣本和SRA文件GEO：http://www.ncbi.nlm.nih.gov/geo.（6個RNA來源），GSE18508（從原文獻中查到GSE號）為csv文件：
> library(limma)> GEO <- readTargets("GEO-samples.csv", sep=",")> GEO GEO Title Pasilla1 GSM461176 S2_DRSC_Untreated-1 Normal2 GSM461177 S2_DRSC_Untreated-3 Normal3 GSM461178 S2_DRSC_Untreated-4 Normal4 GSM461179 S2_DRSC_CG8144_RNAi-1 Down5 GSM461180 S2_DRSC_CG8144_RNAi-3 Down6 GSM461181 S2_DRSC_CG8144_RNAi-4 Down有三種未經處理的生物樣品，其中ps應在正常水平上表達，三種經過處理的生物樣品，其中ps應在表達低水平。
當GEO記錄樣本信息時，序列數據文件實際上保存在NCBI短讀存檔文件(SRA)中。RNA樣本在Illumina基因組分析儀II上進行測序。其中幾個樣本使用了多次測序，共得到20個SRA文件：
> SRA <- readTargets("SRA-Files.csv", sep=",")> SRA SRA GEO Title RunDate FlowCellID Type ReadLength1 SRR031708 GSM461176 S2_DRSC_Untreated-1 7/15/08 308T2AAXX SE 452 SRR031709 GSM461176 S2_DRSC_Untreated-1 7/15/08 308T2AAXX SE 453 SRR031710 GSM461176 S2_DRSC_Untreated-1 8/15/08 30AYWAAXX SE 454 SRR031711 GSM461176 S2_DRSC_Untreated-1 8/15/08 30AYWAAXX SE 455 SRR031712 GSM461176 S2_DRSC_Untreated-1 8/15/08 30AYWAAXX SE 456 SRR031713 GSM461176 S2_DRSC_Untreated-1 8/15/08 30AYWAAXX SE 457 SRR031714 GSM461177 S2_DRSC_Untreated-3 11/14/08 30MNEAAXX PE 378 SRR031715 GSM461177 S2_DRSC_Untreated-3 12/23/08 30M2BAAXX PE 379 SRR031716 GSM461178 S2_DRSC_Untreated-4 12/23/08 30M2BAAXX PE 3710 SRR031717 GSM461178 S2_DRSC_Untreated-4 12/23/08 30M2BAAXX PE 3711 SRR031718 GSM461179 S2_DRSC_CG8144_RNAi-1 7/15/08 308T2AAXX SE 4512 SRR031719 GSM461179 S2_DRSC_CG8144_RNAi-1 7/18/08 308UEAAXX SE 4513 SRR031720 GSM461179 S2_DRSC_CG8144_RNAi-1 8/15/08 30AYWAAXX SE 4514 SRR031721 GSM461179 S2_DRSC_CG8144_RNAi-1 8/15/08 30AYWAAXX SE 4515 SRR031722 GSM461179 S2_DRSC_CG8144_RNAi-1 8/15/08 30AYWAAXX SE 4516 SRR031723 GSM461179 S2_DRSC_CG8144_RNAi-1 8/21/08 308A0AAXX SE 4517 SRR031724 GSM461180 S2_DRSC_CG8144_RNAi-3 12/23/08 30M2BAAXX PE 3718 SRR031725 GSM461180 S2_DRSC_CG8144_RNAi-3 12/23/08 30M2BAAXX PE 3719 SRR031726 GSM461181 S2_DRSC_CG8144_RNAi-4 12/23/08 30M2BAAXX PE 3720 SRR031727 GSM461181 S2_DRSC_CG8144_RNAi-4 12/23/08 30M2BAAXX PE 37##所有在2008年7月和8月進行的測序都使用單端(SE)測序，其鹼基對讀數為45，而後期的測序使用雙端(PE)測序，其鹼基對讀數為37 bp。15.2.3 將reads映射到參照基因組利用Rsubread將reads映射到參考D. melanogaster的基因組：
首先，SRA格式文件需要轉換成FASTQ格式，這是表示高通量測序數據的實際標準。這是使用SRA Toolkit 2.3.4版完成的。具體來說，Unix SRA工具包命令fastq-dump --split-spot --split-files 用於為每個單端SRA文件生成一個FASTQ文件，為每個雙端SRA文件生成兩個FASTQ文件。
第二，代表D.melanogaster的基因組的FASTA文件是從ftp://ftp.ncbi.nlm.nih.gov/genomes/Drosophila_melanogaster/RELEASE_5_48 下載的。染色體臂對應五個FASTA文件：
> DM <- readTargets("FASTA-GFF-files.csv",sep=",")> DM FASTA GFF Chromosome1 CHR_2/NT_033778.fna CHR_2/NT_033778.gff Chr2R2 CHR_2/NT_033779.fna CHR_2/NT_033779.gff Chr2L3 CHR_3/NT_033777.fna CHR_3/NT_033777.gff Chr3R4 CHR_3/NT_037436.fna CHR_3/NT_037436.gff Chr3L5 CHR_4/NC_004353.fna CHR_4/NC_004353.gff Chr46 CHR_X/NC_004354.fna CHR_X/NC_004354.gff ChrX五個*.fna文件被連接到一個單獨的fasta格式文件Drosophila ref.fasta中，並使用Rsubread構建一個索引文件。
> library(Rsubread)> buildindex(basename="index",reference="Drosophila_ref.fasta")第三，reads與基因組對齊。用下面命令映射SE讀數：
> align("index",SE_fastq_files,input_format="FASTQ")其中，SE fastq文件是一個包含SE fastq文件名稱的字符向量。映射PE讀數：
> align("index",PE_fastq_files1,PE_fastq_files2,input_format="FASTQ")其中，PE fastq files1和PE fastq files2是包含PE fastq文件名稱的向量。這些命令生成包含對齊讀數(SE)或片段(PE)的BAM文件。
15.2.4 讀取外顯子計數接下來，我們計算每個基因的注釋外顯子重疊的讀數或片段數。包含基因和外顯子注釋的GFF文件從ftp://ftp.ncbi.nlm.nih. gov/genomes/Drosophila_melanogaster/RELEASE_5_48下載。每個染色體臂都有一個GFF文件，如上所示。
五個*.gff文件連接到一個文件中，並刪除相同外顯子的重複外顯子實例(相同的開始和停止位置)，以創建一個名為unique.gff的開始/停止位置的數據。SE讀數由下面命令進行計數：
> fc_SE <- featureCounts(SE_bam_files, annot.ext="unique.gff",+ isGTFAnnotationFile=TRUE,+ GTF.featureType="exon", GTF.attrType="ID",+ useMetaFeatures=FALSE, allowMultiOverlap=TRUE)SE bam files 是SE讀數名稱的bam文件的向量。PE bam files 是PE讀數名稱的bam文件的向量。PE讀數由下面命令進行計數：
> fc_PE <- featureCounts(PE_bam_files, annot.ext="unique.gff",+ isGTFAnnotationFile=TRUE,+ GTF.featureType="exon", GTF.attrType="ID",+ useMetaFeatures=FALSE, allowMultiOverlap=TRUE,+ isPairedEnd=TRUE)15.2.5 匯總DGEList和計數技術重複的總數#匯總到edgeR DGEList對象> library(edgeR)> y.all <- DGEList(counts=cbind(fc_SE$counts,fc_PE$counts), genes=fc_SE$annotation)> dim(y.all)[1] 74184 20> head(y.all$genes GeneID Chr Start End Strand Length138088 30970 NC_004354.3 138094 139379 - 1286138087 30970 NC_004354.3 139445 139611 - 167138089 30970 NC_004354.3 139445 139889 - 445138086 30970 NC_004354.3 139713 139889 - 177138091 30971 NC_004354.3 140011 141629 + 1619138092 30971 NC_004354.3 142415 144271 + 1857該注釋包括Entrez ID和每個外顯子的長度、染色體和起始、終止位置。回存到原始SRA的順序：
> y.all <- y.all[,SRA$SRA]然後我們通過將計數重複和計數相加把數據摺疊起來，這樣每個GEO樣本就有一個單獨的列：
> y <- sumTechReps(y.all, ID=SRA$GEO)> y$samples group lib.size norm.factorsGSM461176 1 31007529 1GSM461177 1 13040952 1GSM461178 1 15030819 1GSM461179 1 28143539 1GSM461180 1 14901292 1GSM461181 1 16264066 1> colnames(y) <- c("N1","N3","N4","D1","D3","D4")15.2.6 基因注釋果蠅的基因注釋文件：ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Invertebrates。
現在，我們將選中的注釋列添加到DGEList對象中：
``
> ncbi.L1 <- readLines("Drosophila_melanogaster.gene_info", n = 1)> ncbi.colname <- unlist(strsplit(substring(ncbi.L1, 10, 234), 』 』))> ncbi <- read.delim("Drosophila_melanogaster.gene_info",+ skip=1, header=FALSE, stringsAsFactors=FALSE)> colnames(ncbi) <- ncbi.colname> m <- match(y$genes$GeneID, ncbi$GeneID)> y$genes$Chr <- ncbi$chromosome[m]> y$genes$Symbol <- ncbi$Symbol[m]> y$genes$Strand <- NULL> head(y$genes) GeneID Chr Start End Length Symbol138088 30970 X 138094 139379 1286 CG3038138087 30970 X 139445 139611 167 CG3038138089 30970 X 139445 139889 445 CG3038138086 30970 X 139713 139889 177 CG3038138091 30971 X 140011 141629 1619 G9a138092 30971 X 142415 144271 1857 G9a15.2.7 過濾保留至少3個樣本中有超過1個cpm的外顯子。
> isexpr <- rowSums(cpm(y) > 1) >=3> y <- y[isexpr,,keep.lib.sizes=FALSE]15.2.8 Scale 標準化應用TMM標準化
> y <- calcNormFactors(y)> y$sample group lib.size norm.factorsN1 1 30872843 0.955N3 1 12962245 1.031N4 1 14908555 0.976D1 1 27989806 1.005D3 1 14760887 1.022D4 1 16172265 1.014MDS圖顯示了Pasilla down與正常樣本的明顯分離，但也顯示了與測序類型和日期相關的批次效應：
15.2.9 線性模型設置設計矩陣：
> Batch <- factor(c(1,3,4,1,3,4))> Pasilla <- factor(GEO$Pasilla,levels=c("Normal","Down"))> design <- model.matrix(~ Batch + Pasilla)應用voom將讀數轉換為log2-cpm，並賦予相應的權重：
> v <- voom(y,design,plot=TRUE)#線性模型：> fit <- lmFit(v, design)#檢測表達差異的外顯子> fit.de <- eBayes(fit, robust=TRUE)> topTable(fit.de, coef="PasillaDown") GeneID Chr Start End Length Symbol logFC AveExpr t P.Value adj.P.Val B150709 32007 X 10674926 10676128 1203 sesB -3.26 6.41 -32.1 4.45e-14 8.35e-10 22.3150713 32007 X 10675026 10676128 1103 sesB -3.26 6.41 -32.1 4.48e-14 8.35e-10 22.396435 43230 3R 22694960 22695852 893 BM-40-SPARC -2.47 10.50 -28.2 2.51e-13 2.61e-09 20.9150697 32008 X 10672987 10673728 742 Ant2 2.85 5.53 28.0 2.80e-13 2.61e-09 20.696434 43230 3R 22695915 22696094 180 BM-40-SPARC -2.27 8.11 -26.9 4.72e-13 3.52e-09 20.396433 43230 3R 22697648 22697717 70 BM-40-SPARC -2.17 6.27 -23.5 2.70e-12 1.46e-08 18.6107856 44030 3L 2561932 2562843 912 msn -2.46 5.12 -23.5 2.80e-12 1.46e-08 18.570750 44258 3R 5271691 5272628 938 ps -2.27 5.53 -23.1 3.38e-12 1.46e-08 18.411333 44548 2R 6407125 6408782 1658 lola 2.25 5.73 23.1 3.54e-12 1.46e-08 18.3150702 32008 X 10674230 10674694 465 Ant2 2.97 3.89 21.9 6.82e-12 2.34e-08 17.4> summary(decideTests(fit.de)) (Intercept) Batch3 Batch4 PasillaDownDown 443 5380 6070 2062NotSig 3388 26386 25184 33346Up 33427 5492 6004 185015.2.10 可變剪接Pasilla down與normal之間外顯子保留的差異檢測：
> ex <- diffSplice(fit[,"PasillaDown"], geneid = "GeneID", exonid = "Start")Total number of exons: 37258Total number of genes: 8192Number of genes with 1 exon: 1619Mean number of exons in a gene: 5Max number of exons in a gene: 56顯示差異剪接的基因：
> topSplice(ex, test="simes", n=20) GeneID Chr Symbol NExons P.Value FDR11214 44548 2R lola 30 5.52e-32 3.63e-2895956 44448 3R scrib 35 1.99e-20 4.41e-17141235 45320 X trol 44 2.01e-20 4.41e-17
> topSplice(ex, test="F", n=20) GeneID Chr Symbol NExons F P.Value FDR141235 45320 X trol 44 30.87 1.93e-40 1.27e-3611214 44548 2R lola 30 41.80 4.81e-33 1.58e-2995956 44448 3R scrib 35 16.20 1.53e-23 3.35e-20##由於篇幅限制，這裡只顯示了top3，因為與下面的三幅圖像相關。我們可以通過繪製所有的外顯子顯示差異最大的剪接基因：
> plotSplice(ex, geneid="lola", genecol="Symbol")> plotSplice(ex, geneid="scrib", genecol="Symbol")> plotSplice(ex, geneid="trol", genecol="Symbol")我們可以看到，當Pasilla基因表達下降時，trol基因右端一個由5或6個外顯子組成的基因塊相對丟失，代價是較早的一個由12個外顯子組成的基因塊在Pasilla基因表達下降時相對高表達。基因前半部分的大多數外顯子的行為相似。