RNA-seq測序數(shù)據(jù)分析

RNA測序數(shù)據(jù)分析RNAdataanalysisRNAbigdataanalysisRNAsequencing:Transcriptome目錄

RNA-seq簡述

RNA-seq常規(guī)分析

RNA-seq高級分析

RNA-seq拓展RNA-seq簡述RNA-seq簡述轉(zhuǎn)錄組學(xué)The

transcriptome

isthesetofall

RNA

moleculesinonecellorapopulationof

cells.Itissometimesusedtoreferto

allRNAs,orjust

mRNA,dependingontheparticularexperiment.WangZ.NatRevGenet,2009.RNA-seq簡述轉(zhuǎn)錄組學(xué)WangZ.NatRevGenet,2009.RNA-seq簡述轉(zhuǎn)錄組學(xué)RNA-seq簡述文庫制備RNA-seq簡述測序流程MARTINJA.NatRevGenet,2011.RNA-seq簡述術(shù)語RNA-seq簡述數(shù)據(jù)分析目錄RNA-seq簡述

RNA-seq常規(guī)分析

RNA-seq高級分析

RNA-seq拓展RNA-seq常規(guī)分析

RNA-seq前期處理

基因定量和表達差異性分析

RNA-seq組裝RNA-seq常規(guī)分析前期處理RNA-seq常規(guī)分析前期處理RNA-seq常規(guī)分析前期處理RNA-seq常規(guī)分析前期處理RNA-seq常規(guī)分析前期處理下載測序數(shù)據(jù)以NCBI

SRA數(shù)據(jù)庫為例：下載測序數(shù)據(jù)下載測序數(shù)據(jù)SRAtoolkit工具解壓原始.sra文件：軟件下載：/Traces/sra/sra.cgi?view=softwarewget/sra/sdk/current/sratoolkit.current-centos_linux64.tar.gz下載測序數(shù)據(jù)SRAtoolkit安裝：tar

xzf

sratoolkit.current-centos_linux64.tar.gz解壓直接使用即可，里面有一大堆的軟件，針對不同的測序儀，不同的數(shù)據(jù)下載測序數(shù)據(jù)SRAtoolkit使用：使用fastq-dump命令將SRA格式轉(zhuǎn)換成FASTQ格式：基本用法：fastq-dumpSRR123456.sra參數(shù)說明：-O指定輸出目錄–split-3對于雙末端測序，一定使用這個參數(shù)，將一個*.sra文件轉(zhuǎn)換成*_1.fastq和*_2.fastq兩個文件，否則就會轉(zhuǎn)換成一個*.fastq文件，后續(xù)的分析會出錯。下載測序數(shù)據(jù)SRAtoolkit使用：使用fastq-dump命令將SRA格式轉(zhuǎn)換成FASTQ格式（批量）：foriin*sradoecho$i~

/fastq-dump--split-3$idone測序數(shù)據(jù)預(yù)處理為什么要對reads進行前期處理？1，測序儀器硬件原因產(chǎn)生的信號強度極端的reads；2，測序熒光強度不夠造成reads中含有的模糊的N堿基；3，reads中含有接頭序列（adapter）；4，總體質(zhì)量偏低的reads。即Q>=20堿基比例小于50%的

reads，其中Q=-10logerror_ratio測序數(shù)據(jù)預(yù)處理RawData去除reads中堿基為N的堿基數(shù)達到一定比例的reads去除低質(zhì)量堿基數(shù)目達到一定程度的reads去除adapter污染的reads去除duplicationreads(差異性為0)UsedData測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量FastQC工具：KatherineT.A.PlosOne.2012.下載安裝：可以從FastQC官網(wǎng)上下載，軟件http://www.bioinformatics.babraham.ac.uk/projects/fastqc/1，軟件支持多個系統(tǒng)平臺，包括windows版、Linux和MacOS。注意：windows版本需要Java運行環(huán)境，處理大數(shù)據(jù)也比較消耗內(nèi)存；2，Linux版本。使用unzip命令來進行解壓縮。然后進入解壓縮文件，fastqc文件即是主程序。我們使用chmodu+x命令修改為可執(zhí)行權(quán)限。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量KatherineT.A.PlosOne.2012.fastqc[-ooutputdir][--(no)extract][-ffastq|bam|sam][-ccontaminantfile]seqfile1..seqfileN-o用來指定輸出文件的所在目錄，注意是不能自動新建目錄的。輸出的結(jié)果是.zip文件，默認自動解壓縮，命令里加上—(no)extract則不解壓縮。-f用來強制指定輸入文件格式，默認會自動檢測。-c用來指定一個contaminant文件，fastqc會把overrepresentedsequences往這個contaminant文件里搜索。contaminant文件的格式是"Name\tSequences”。加上-q會進入沉默模式，即不出現(xiàn)下面的提示：?Startedanalysisoftarget.fqApprox5%completefor

target.fq?Approx10%completefor

target.fq如果輸入的fastq文件名是target.fq，fastqc的輸出的壓縮文件將是target.fq_fastqc.zip。解壓后，查看html格式的結(jié)果報告。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量結(jié)果分為綠色的"PASS"，黃色的"WARN"和紅色的"FAIL"。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量1Basicstatistics測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量2Perbasesequencequalityquality就是Fred值，-10*log10(p)，p為測錯的概率。所以一條reads某位置出錯概率為0.01時，其quality就是20。橫軸代表位置，縱軸quality。紅色表示中位數(shù)，黃色是25%-75%區(qū)間，觸須是10%-90%區(qū)間，藍線是平均數(shù)。若任一位置的下四分位數(shù)低于10或中位數(shù)低于25，報“WARN”；若任一位置的下四分位數(shù)低于5或中位數(shù)低于20，報“FAIL”。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量3PerSequenceQualityScores橫軸為quality，縱軸是reads數(shù)目。當(dāng)峰值小于27（錯誤率0.2%）時報"WARN"，當(dāng)峰值小于20（錯誤率1%）時報"FAIL"。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量4PerBaseSequenceContent橫軸為位置，縱軸為百分比。正常情況下四種堿基的出現(xiàn)頻率應(yīng)該是接近的，而且沒有位置差異。因此好的樣本中四條線應(yīng)該平行且接近。當(dāng)部分位置堿基的比例出現(xiàn)bias時，即四條線在某些位置紛亂交織，往往提示我們有overrepresentedsequence的污染。當(dāng)所有位置的堿基比例一致的表現(xiàn)出bias時，即四條線平行但分開，往往代表文庫有bias（建庫過程或本身特點），或者是測序中的系統(tǒng)誤差。當(dāng)任一位置的A/T比例與G/C比例相差超過10%，報"WARN"；當(dāng)任一位置的A/T比例與G/C比例相差超過20%，報"FAIL"。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量5PerBaseGCContent如果建庫足夠均勻，reads的每個位置應(yīng)當(dāng)是沒有差異的，所以GC含量的線應(yīng)當(dāng)平行于X軸，反映樣品（基因組、轉(zhuǎn)錄組等）的GC含量。當(dāng)部分位置GC含量出現(xiàn)bias時，往往提示我們有overrepresentedsequence的污染。當(dāng)所有位置的GC含量一致的表現(xiàn)出bias時，往往代表文庫有bias(建庫過程或本身特點)，或者是測序中的系統(tǒng)誤差。當(dāng)任一位置的GC含量偏離均值的5%時，報"WARN"；當(dāng)任一位置的GC含量偏離均值的10%時，報"FAIL"。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量6PerSequenceGCContent紅線是實際情況，藍線是理論分布（正態(tài)分布，均值不一定在50%，而是由平均GC含量推斷的）。曲線形狀的偏差往往是由于文庫的污染或是部分reads構(gòu)成的子集有偏差（overrepresentedreads）。形狀接近正態(tài)但偏離理論分布的情況提示我們可能有系統(tǒng)偏差。偏離理論分布的reads超過15%時，報"WARN"；偏離理論分布的reads超過30%時，報"FAIL"。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量7PerBaseNContent正常情況下N的比例是很小的，所以圖上常?？吹揭粭l直線，但放大Y軸之后會發(fā)現(xiàn)還是有N的存在，這不算問題。當(dāng)Y軸在0%-100%的范圍內(nèi)也能看到“鼓包”時，說明測序系統(tǒng)出了問題。當(dāng)任意位置的N的比例超過5%，報"WARN"；當(dāng)任意位置的N的比例超過20%，報"FAIL"。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量8SequenceLengthDistribution當(dāng)reads長度不一致時報"WARN"；當(dāng)有長度為0的read時報“FAIL”。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量9DuplicateSequences統(tǒng)計序列完全一樣的reads的頻率。測序深度越高，越容易產(chǎn)生一定程度的duplication，這是正常的現(xiàn)象，但如果duplication的程度很高，就提示我們可能有bias的存在（如建庫過程中的PCRduplication）。橫坐標是duplication的次數(shù)，縱坐標是duplicatedreads的數(shù)目，以uniquereads的總數(shù)作為100%。fastqc中用fq數(shù)據(jù)的前200,000條reads統(tǒng)計其在全部數(shù)據(jù)中的重復(fù)情況。重復(fù)數(shù)目大于等于10的reads被合并統(tǒng)計，這也是為什么我們看到上圖的最右側(cè)略有上揚。當(dāng)非unique的reads占總數(shù)的比例大于20%時，報“WARN”；當(dāng)非unique的reads占總數(shù)的比例大于50%時，報“FAIL”。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量10OverrepresentedSequences如果有某個序列大量出現(xiàn)，就叫做over-represented。fastqc的標準是占全部reads的0.1%以上。和上面的duplicateanalysis一樣，為了計算方便，只取了fq數(shù)據(jù)的前200,000條reads進行統(tǒng)計，所以有可能over-representedreads不在里面。而且大于75bp的reads也是只取50bp。如果命令行中加入了-ccontaminantfile，出現(xiàn)的over-representedsequence會從contaminant_file里面找匹配的hit（至少20bp且最多一個mismatch），可以給我們一些線索。當(dāng)發(fā)現(xiàn)超過總reads數(shù)0.1%的reads時報”WARN“，當(dāng)發(fā)現(xiàn)超過總reads數(shù)1%的reads時報”FAIL“。測序數(shù)據(jù)預(yù)處理檢查序列質(zhì)量11OverrepresentedKmers如果某k個bp的短序列在reads中大量出現(xiàn)，其頻率高于統(tǒng)計期望的話，fastqc將其記為over-representedk-mer。默認的k=5，可以用-k--kmers選項來調(diào)節(jié)，范圍是2-10。出現(xiàn)頻率總體上3倍于期望或是在某位置上5倍于期望的k-mer被認為是over-represented。fastqc除了列出所有over-representedk-mers，還會把前6個的perbasedistribution畫出來。當(dāng)有出現(xiàn)頻率總體上3倍于期望或是在某位置上5倍于期望的k-mer時，報”WARN“；當(dāng)有出現(xiàn)頻率在某位置上10倍于期望的k-mer時報"FAIL"。測序數(shù)據(jù)預(yù)處理去除低質(zhì)量序列Trimmomatic工具：Trim

readsBolger.A.M.Bioinformatics.2014.Trimmomatic使用JAVA運行，速度快下載地址:/cms/index.php?page=trimmomatic下載binary文件，直接執(zhí)行命令即可測序數(shù)據(jù)預(yù)處理去除低質(zhì)量序列Trimmomatic工具：Trim

readsBolger.A.M.Bioinformatics.2014.測序數(shù)據(jù)預(yù)處理去除低質(zhì)量序列Bolger.A.M.Bioinformatics.2014.PairedEnd:java-jartrimmomatic-0.35.jarPE–thread8-phred33input_forward.fq.gzinput_reverse.fq.gzoutput_forward_paired.fq.gzoutput_forward_unpaired.fq.gzoutput_reverse_paired.fq.gzoutput_reverse_unpaired.fq.gzILLUMINACLIP:TruSeq3-PE.fa:2:30:10LEADING:3TRAILING:3SLIDINGWINDOW:4:15MINLEN:36PE/SE設(shè)定對Paired-End或Single-End的reads進行處理，其輸入和輸出參數(shù)稍有不一樣。-threads設(shè)置多線程運行數(shù)-phred33設(shè)置堿基的質(zhì)量格式，可選pred64ILLUMINACLIP:TruSeq3-PE.fa:2:30:10切除adapter序列。參數(shù)后面分別接adapter序列的fasta文件：允許的最大mismatch數(shù)：palindrome模式下匹配堿基數(shù)閾值：simple模式下的匹配堿基數(shù)閾值。LEADING:3切除首端堿基質(zhì)量小于3的堿基TRAILING:3切除尾端堿基質(zhì)量小于3的堿基SLIDINGWINDOW:4:15

Performaslidingwindowtrimming。Windows的size是4個堿基，其平均堿基質(zhì)量小于15，則切除。MINLEN:50

最小的reads長度CROP:<length>

保留reads到指定的長度HEADCROP:<length>

在reads的首端切除指定的長度TOPHRED33將堿基質(zhì)量轉(zhuǎn)換為pred33格式TOPHRED64將堿基質(zhì)量轉(zhuǎn)換為pred64格式測序數(shù)據(jù)預(yù)處理去除低質(zhì)量序列Bolger.A.M.Bioinformatics.2014.SingleEnd:java-jartrimmomatic-0.35.jarSE-phred33input.fq.gzoutput.fq.gzILLUMINACLIP:TruSeq3-SE:2:30:10LEADING:3TRAILING:3SLIDINGWINDOW:4:15MINLEN:36Thiswillperformthefollowing:Removeadapters(ILLUMINACLIP:TruSeq3-PE.fa:2:30:10)RemoveleadinglowqualityorNbases(belowquality3)(LEADING:3)RemovetrailinglowqualityorNbases(belowquality3)(TRAILING:3)Scanthereadwitha4-basewideslidingwindow,cuttingwhentheaveragequalityperbasedropsbelow15(SLIDINGWINDOW:4:15)Dropreadsbelowthe36baseslong(MINLEN:36)測序數(shù)據(jù)預(yù)處理去除低質(zhì)量序列去除低質(zhì)量reads一般經(jīng)驗：前20bp沒有N質(zhì)量值低于10的堿基比例不超過10%質(zhì)量值低于13的堿基比例不超過20%reads平均質(zhì)量值大于等于20RNA-seq常規(guī)分析RNA-seq前期處理

基因定量和表達差異性分析

RNA-seq組裝RNA-seq常規(guī)分析組裝轉(zhuǎn)錄組組裝：1，基于參考基因組組裝2，從頭組裝轉(zhuǎn)錄組3，二者相結(jié)合RNA-seq常規(guī)分析組裝基于參考基因組MARTINJA.NatRevGenet,2011.RNA-seq常規(guī)分析組裝基于參考基因組MARTINJA.NatRevGenet,2011.ScriptureCufflinksRNA-seq常規(guī)分析組裝基于參考基因組優(yōu)點一個大的拼接問題（上百萬測序序列）轉(zhuǎn)化小的拼接問題拼接敏感性高且能組裝出低豐度的轉(zhuǎn)錄本不足太依賴于參考基因組的質(zhì)量橫跨內(nèi)含子的spliced測序序列容易丟失比對工具的隨機比對（randomassignment）

MARTINJA.NatRevGenet,2011.RNA-seq常規(guī)分析組裝MARTINJA.NatRevGenet,2011.從頭組裝RNA-seq常規(guī)分析組裝MARTINJA.NatRevGenet,2011.Trintiy

SOAPdenovo-TransTrans-ABySSVelvet-OasesIDBA-TranRNA-seq常規(guī)分析組裝MARTINJA.NatRevGenet,2011.基于從頭組裝優(yōu)點不依賴于參考基因組不需要一些比對工具可識別新可變剪切點不足在識別高等真核生物可變剪切同源異構(gòu)體存在挑戰(zhàn)全長轉(zhuǎn)錄組組裝過程需要高的序列測序深度對測序錯誤敏感RNA-seq常規(guī)分析組裝MARTINJA.NatRevGenet,2011.Reference-baseddenovo二者組合組裝RNA-seq常規(guī)分析組裝RNA-seq常規(guī)分析組裝MARTINJA.NatRevGenet,2011.RNA-seq常規(guī)分析組裝MARTINJA.NatRevGenet,2011.RNA-seq常規(guī)分析組裝長度分布N50N90靈敏度特異度ROCRNA-seq常規(guī)分析組裝TranscriptassemblyperformanceSteijgerT.NatMethods,2013.RNA-seq常規(guī)分析組裝評估TranscriptassemblyperformanceSteijgerT.NatMethods,2013.RNA-seq常規(guī)分析組裝RNA-seq常規(guī)分析組裝序列比對BWAHengL.Bioinformatics.2009.下載地址：http://安裝：#installingBWAtarjxvf~/software/bwa-0.7.10.tar.bz2-C$yourPATHmakeecho'PATH=$PATH$yourPATH'>>~/.bashrcsource~/.bashrc序列比對BWAHengL.Bioinformatics.2009.BWA使用步驟：使用BWA構(gòu)建參考基因組的index數(shù)據(jù)庫BWA-MEM比對（BWA-SW和BWA-backtrack用的少）BWA-backtrack,BWA-SWandBWA-MEM.ThefirstalgorithmisdesignedforIlluminasequencereadsupto100bp,whiletheresttwoforlongersequencesrangedfrom70bpto1Mbp.BWA-MEMandBWA-SWsharesimilarfeaturessuchaslong-readsupportandsplitalignment,butBWA-MEM,whichisthelatest,isgenerallyrecommendedforhigh-qualityqueriesasitisfasterandmoreaccurate.BWA-MEMalsohasbetterperformancethanBWA-backtrackfor70-100bpIlluminareads.序列比對BWAHengL.Bioinformatics.2009.1,使用BWA構(gòu)建參考基因組的index數(shù)據(jù)庫序列比對BWAHengL.Bioinformatics.2009.2,尋找SA

coordinates如果是pair-end數(shù)據(jù)（leftRead.fastq和rightRead.fastq）兩個文件分別處理bwaalnreference.faleftRead.fastq>leftRead.saibwaalnreference.farightRead.fastq>rightRead.sai若是single-end數(shù)據(jù)，直接單獨處理bwaalnreference.fasingleRead.fastq>singleRead.sai序列比對BWAHengL.Bioinformatics.2009.3,轉(zhuǎn)換SAcoordinates輸出為sam若是paired-end:bwasampe-fpair-end.samreference.faleftRead.sairightRead.saileftRead.fastqrightread.fastq若是single-end:bwasamse-fsingle.samreference.fasingle.saisingle.fastq序列比對BWAHengL.Bioinformatics.2009.可選擇的：bwamemref.fareads.fq>aln-se.sambwamemref.faread1.fqread2.fq>aln-pe.sambwabwaswref.falong_read.fq>aln.samRNA-seq常規(guī)分析RNA-seq組裝

基因定量和表達差異性分析RNA-seq前期處理RNA-seq常規(guī)分析定量差異性分析RNA-seq常規(guī)分析定量差異性分析RNA-seq常規(guī)分析定量差異性分析RPKM(ReadsPerKilobaseMillion),FPKM(FragmentsPerKilobaseMillion)和TPM(TransciptsPerMillion)作為標準化數(shù)值。RNA-seq常規(guī)分析定量差異性分析RPKM(假設(shè)以10為單位)RNA-seq常規(guī)分析定量差異性分析TPM(假設(shè)以10為單位)RNA-seq常規(guī)分析定量差異性分析TPMRPKM每個樣本的TPM的總和是相同的，這就意味著TPM數(shù)值能體現(xiàn)出比對上某個基因的reads的比例，使得該數(shù)值可以直接進行樣本間的比較。RNA-seq常規(guī)分析定量差異性分析RNA-seq常規(guī)分析定量差異性分析RNA-seq常規(guī)分析定量差異性分析RNA-seq常規(guī)分析定量差異性分析RNA-seq常規(guī)分析定量差異性分析MihaelaP.NatProtocal,2016.RNA-seq常規(guī)分析定量差異性分析DaehwanK.NatMethods,2015.HisatRNA-seq常規(guī)分析定量差異性分析MihaelaP.NatBiotechnology,2015.StringTie轉(zhuǎn)錄組數(shù)據(jù)比對BowtieKatherineT.A.PlosOne.2012.Bowtie下載：/index.shtmlunzipbowtie-0.12.7-macos-10.5-x86_64.zipcdbowtie-0.12.7

1對Reference文件(GENOME.fa)建庫bowtie-buildGENOME.faGENOME (網(wǎng)站上可以直接下載已indexed的庫文件)2使用：

bowtie[options]*<ebwt>{-1<m1>-2<m2>|--12<r>|<s>}[<hit>]

ebwt:indexfile-1-2:paired-end/matepairedreads--12:unpairedandpaired-endreadshit:outputmappedfile轉(zhuǎn)錄組數(shù)據(jù)比對TopHatKatherineT.A.PlosOne.2012.TopHat(/)

$tarzxvftophat-1.3.2.OSX_x86_64.tar.gz

$cdtophat-1.3.2.OSX_x86_64在使用TopHat前，必須將Bowtie的可執(zhí)行文件的目錄輸出到PATH變量中去exportPATH=$PATH:/share/sbin/bowtie 確保TopHat可以運行bowtie,bowtie-inspect以及bowtie-build。還需要下載安裝samtools。TopHat的使用范例：tophat[options]*<ebwt_base><reads1_1[,...,readsN_1]>[reads1_2,...readsN_2]

轉(zhuǎn)錄組組裝CufflinksCufflinks安裝(http://cuff/)

tarzxvfcuffinks-1.2.1.OSX_x86_64.tar.gz

cdcuffinks-1.2.1.OSX_x86_64

cufflinks套裝有很多，我們主要使用的只有三個cufflinks是用來處理tophat的輸出的bam文件然后輸出gtf文件cuffmerge把多個樣本的gtf文件合并的，主要是測多個樣本可能會需要cuffdiff算出分組的bam文件里面的差異基因。轉(zhuǎn)錄組組裝Cufflinks轉(zhuǎn)錄組組裝轉(zhuǎn)錄組組裝TrinityTrinity()

表達差異分析Cufflinks+CummeRbund

表達差異分析Trinity+edgeR目錄RNA-seq簡述

RNA-seq常規(guī)分析

RNA-seq高級分析

RNA-seq拓展RNA-seq高級分析

可視化Co-expression

分析Gopathway分析RNA-seq高級分析可視化RNA-seq高級分析可視化RNA-seq高級分析可視化RNA-seq高級分析可視化RNA-seq高級分析可視化RNA-seq高級分析可視化散點圖折線圖RNA-seq高級分析可視化箱線圖RNA-seq高級分析可視化RNA-seq高級分析可視化log2(foldchange)RNA-seq高級分析可視化RNA-seq高級分析可視化RNA-seq高級分析

可視化Co-expression

分析

Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析Gopathway分析RNA-seq高級分析

可視化

Co-expression