2017年轉(zhuǎn)錄組入門(mén)合輯

【PPT】文獻(xiàn)報(bào)告-AsurveyofbestpracticesforRNA-seq 講解............................................................................................................................................... 內(nèi) 轉(zhuǎn)錄組相關(guān)的安 講解.............................................................................................................................................準(zhǔn) 其他安裝筆 讀文章拿到數(shù) 講解.............................................................................................................................................文 數(shù)據(jù).............................................................................................................................................其他學(xué)習(xí)筆 了解fastq數(shù) 講解.............................................................................................................................................sra文件轉(zhuǎn)換為fastq格 數(shù)據(jù)校驗(yàn)及備份存 質(zhì)控 了解參考組及注 講解.............................................................................................................................................任務(wù)列 在UCSChg19參考 從gencode數(shù)據(jù)庫(kù)注釋文件，并且用IGV去查看感的的結(jié) 截圖幾個(gè)的IGV可視化結(jié) ENSEMBL，NCBI的gtf 序列比 任務(wù)列 比對(duì).............................................................................................................................................hisat2的用 index文件 比對(duì)、排序、索 質(zhì)量控 載入IGV，截圖幾個(gè).................................................................................................................reads計(jì) 任 代碼記 參考資 個(gè)人筆記平 博 .........................................................................................................................................【PPT】文獻(xiàn)報(bào)告-AsurveyofbestpracticesforRNA-seqdataAsurveyofbestpracticesforRNA-seqdataysis我把它叫做RNA-seq數(shù)據(jù)分析指南。這篇文章是由佛羅里達(dá)大學(xué)等單位的研究在1月26日在GenomeBiology上的，該期刊的影響因子有10.8分。這是這篇文章的通訊作者，應(yīng)該挺靠譜的。新一代技術(shù)在式發(fā)展的同時(shí)，也衍生出許多其他技術(shù)創(chuàng)新。RNA-Seq就是其中之seq可以獲得相當(dāng)驚人的數(shù)據(jù)量，而這恰恰是一柄雙刃劍。豐富的數(shù)據(jù)量蘊(yùn)含著大量的寶貴正因如此，數(shù)據(jù)分析可以說(shuō)是RNA-seq的重中之重。RNA-seq有非常廣泛的應(yīng)用，但沒(méi)有分析策略?，F(xiàn)在人們已經(jīng)了大量的RNA-seq和數(shù)據(jù)分析方案，對(duì)于剛?cè)腴T(mén)的新手來(lái)說(shuō)難免有些無(wú)所適從。這篇文章概述了RNA-seq生物信息學(xué)分析的準(zhǔn)和現(xiàn)有資源，為人們提供了一份RNA-seq數(shù)據(jù)分析指南，可以作為開(kāi)展RNA-seq研究的寶貴參考資料。這份指南覆蓋了RNA-seq數(shù)據(jù)分析的所有主要步驟，比如質(zhì)量控制、讀段比對(duì)、和轉(zhuǎn)錄本定量、差異性表達(dá)、功能分析、融合檢測(cè)、eQTL圖譜分析等等。研究繪制的RNA-seq分析通用路線(xiàn)圖(標(biāo)準(zhǔn)Illumina)，將主要分析步驟分為前期分析、分析從這，橫軸是年份，縱軸是儀的通量，圈里面的數(shù)字代表讀長(zhǎng)?？梢钥吹絻x的通量和讀長(zhǎng)都在增加了，其中PacficBiosciences的三代儀讀長(zhǎng)最長(zhǎng)，可以達(dá)到14K，illumina的儀通量最大，HiSeqXTen的通量可以達(dá)到1.8T?，F(xiàn)在Illumina生從左到右，儀的通量逐漸增大，它們適合不同的樣品和目的。Miseq通量比較低，適合宏組等微生物；Hiseq通量太高了，如果你送去公司，他們一般要20~30這幅圖的橫軸是年份，縱軸是高通量技術(shù)應(yīng)用的代表性文章的量。不同的應(yīng)用技術(shù)用顏色進(jìn)行分類(lèi)，數(shù)據(jù)點(diǎn)的大小跟率（率/月）成正比?？梢钥闯鯮NA-Seq技術(shù)的mRNA在生物內(nèi)RNA的組分中只占很小的一部分，rRNA占絕大多數(shù)。一般說(shuō)RNA-seq指的都是mRNA-seq，后面的流程也都是主要針對(duì)mRNA-seq數(shù)據(jù)分析的。在科學(xué)家們的努力下，可以把那些非編碼RNA提取出來(lái)建庫(kù)，進(jìn)序。對(duì)于Illumina，片段一般小于500bp。確定合適長(zhǎng)度的片段是后續(xù)和分析然后呢，設(shè)定生物學(xué)重復(fù)對(duì)差異的檢出率（真陽(yáng)性率，TPR）的提高具有明顯效果。上數(shù)據(jù)量分配到的生物學(xué)重復(fù)樣本中，差異分析結(jié)果的可靠性在不斷提升。對(duì)于RNA-seq，設(shè)置一些重復(fù)，差異不穩(wěn)定的話(huà)有時(shí)候設(shè)置10個(gè)/20個(gè)都不夠。具體問(wèn)題具體分析??！深度（Sequencingdepth），也叫乘數(shù)，指每個(gè)堿基被的平均次數(shù)，是用來(lái)衡量測(cè)序量的首要參數(shù)。研究表明，增加深度，量從1.6M條reads增加到20M條reads,（75bp）10Mreads時(shí)就已經(jīng)達(dá)到平衡了，80%的雞轉(zhuǎn)錄本被檢測(cè)到。在此基礎(chǔ)上增直觀(guān)一些說(shuō)，如果某個(gè)在RNA-seq結(jié)果顯示差異表達(dá)，但QPCR結(jié)果表明這個(gè)表達(dá)差異不顯著，可以認(rèn)為這個(gè)RNA-seq結(jié)果為假陽(yáng)性；反之，這個(gè)結(jié)果就是真陽(yáng)性。如表2、3所示，在一定的生物學(xué)重復(fù)數(shù)（n）的情況下，隨著單樣本量（Depth）的提高（25%→100%），真陽(yáng)性率（TPR）都只有有限的提高。例如在n=3的情況下，單個(gè)樣本的量從25%提高到100%，TPR僅僅從6.24%提高到8.95%。在表3中，如果Depth等于25%不變，當(dāng)n從2提高到12，TPR的提高則是非常明顯的。因此深度對(duì)結(jié)果改為n*1/n=1,保持不變）。如圖A,灰色實(shí)線(xiàn)代表不同的生物學(xué)重復(fù)數(shù)（n）和單樣本數(shù)據(jù)量高。例如n=2,TPR約為3%，如果n=6,TPR則提高到22%。但單個(gè)樣本的數(shù)據(jù)量不斷降低，TPR的降低十分緩慢。例如，n=3,單個(gè)樣本的數(shù)據(jù)量從15%，TPR9%5%和單樣本量的組合，對(duì)假陽(yáng)性率（FPR）的影響卻較小。如圖B，灰色實(shí)線(xiàn)代表不同生物學(xué)重復(fù)數(shù)（n）和單樣本數(shù)據(jù)量（1/n）組合的情況下，真陽(yáng)性率（FPR）的變化。雖然n從2變化到96，F(xiàn)PR基本沒(méi)有太大變化。從圖中很容易發(fā)現(xiàn)，基于負(fù)二項(xiàng)分布的差異分析檢驗(yàn)（Pvalue），F(xiàn)PR對(duì)生物學(xué)重復(fù)數(shù)和單個(gè)樣本數(shù)據(jù)量均不敏感，始終保持低于0.1%水平?；蛘哒f(shuō)，這個(gè)算法對(duì)FPR的控制還隨著單價(jià)的下降，目前市場(chǎng)上RNA-seq類(lèi)項(xiàng)目的單樣本量正在不斷提高。以2G，PE100的表達(dá)譜項(xiàng)目為例，其對(duì)應(yīng)的量為20M條reads。如果一條長(zhǎng)度為1kbp的低表達(dá)的表達(dá)量為RPKM=0.5，其理論上可以檢測(cè)到的reads數(shù)為20×0.5=10。所以低豐度的檢測(cè)，對(duì)RNA-seq這個(gè)技術(shù)來(lái)說(shuō)并非最大問(wèn)題。的表達(dá)量變化程度，可以使用Qpcr來(lái)驗(yàn)證。但往往也使用所有差異來(lái)統(tǒng)計(jì)某些規(guī)律。例如使用差異的pathway富集分析來(lái)尋找與性狀相關(guān)的pathway。如果在全局水平的差異集并不可靠，那么pathway富集分析得出的結(jié)論的可靠性自然也受到影響。而全局水平的差異數(shù)量巨大，是難以使用Qpcr驗(yàn)證的。因此，定量以及差異分析的準(zhǔn)確性是在RNA-seq中更值得關(guān)心的問(wèn)題。run跑不完，為了避免技術(shù)誤差造成太大的實(shí)驗(yàn)誤差，要把樣品隨機(jī)分配到每個(gè)批次或runs中;(到底怎么設(shè)計(jì)，要一下！！)如果你的樣品是多樣品混合，每個(gè)樣品要單獨(dú)加上，每個(gè)lanes要保證足夠的測(cè)序深度，為了保證所有的樣品在每個(gè)lane中都有。如果送給公司去做的話(huà)，要選擇建庫(kù)llumina平臺(tái)用FastQC看；其他平臺(tái)的數(shù)據(jù)用NGSQC。一般會(huì)有原始數(shù)據(jù)的序列質(zhì)量，GC含量，存在的接頭以及K-mers子串圖并且重復(fù)序列太多的reads。在下降，因此數(shù)據(jù)質(zhì)量逐漸降低乃是自然趨勢(shì)。常用的數(shù)據(jù)過(guò)濾的有FASTX-Reads比對(duì)后的質(zhì)量控制(評(píng)估比對(duì)質(zhì)量的指標(biāo))：比對(duì)上的reads占總reads的百分比；Reads比對(duì)到外顯子和參考鏈上的覆蓋度是否一致；比對(duì)到組序列：多重比對(duì)reads？錄本被定量以后，應(yīng)該看一下GC含量和長(zhǎng)度偏差，確定定量的方法是否適用。把所有樣本的reads混合用于轉(zhuǎn)錄本的拼接。二代的轉(zhuǎn)錄組reads用于拼接還是存在一些問(wèn)題的，最終拼接結(jié)果不太理想。一個(gè)轉(zhuǎn)錄本的拼接結(jié)果會(huì)是10~100contigs。三代是一樣的，F(xiàn)PKM可以轉(zhuǎn)換成TPM。Cufflinks(支持雙端數(shù)據(jù)，并且需要GTF格式的注RNA-seqRNA18~34個(gè)堿基，包含了miRNAs,short-interferingRNAssiRNAs)，PIWI-interactingRNAs(piRNAs)以及其他種類(lèi)的**。sRNA-seqlibrariesarerarelysequencedasdeeplyasregularRNA-seqlibrariesbecauseofalackofcomplexity,withatypicalrangeof2–10millionreads.miRTools2.0atoolforpredictionandprofilingofsRNAspeciesusesbydefaultreadsthatare18–30baseslong5.比對(duì)到參考組上，比對(duì)有：Bowtie2,STAR,orBurrows-WheelerAligner(BWA)PatMaNandMicroRazerSmapshortMoreover,thecombinationofRNA-seqandre-sequencingcanbeusedbothtoremovefalsepositiveswheninferringfusiongenesandto yzecopynumberalterations.Thestatisticallysignificantcorrelationsthatwereobserved,however,accountedforrelativelysmalleffects.(DNAmethylation)一些：CORNA,MMIA,,MAGIA,and代謝組和轉(zhuǎn)錄組數(shù)據(jù)結(jié)合進(jìn)行通路分析，有一些：MassTRIX,Paintomics,VANTEDv2,andSteinerNetRNA-seq技術(shù)已經(jīng)成為轉(zhuǎn)錄組分析的標(biāo)準(zhǔn)方法。其相對(duì)應(yīng)的技術(shù)和數(shù)據(jù)分析工具還在不斷地發(fā)展。對(duì)低表達(dá)的的定量仍是一個(gè)等待解決的問(wèn)題;三代技術(shù)，Smart-seq和Smart-seq2應(yīng)用于轉(zhuǎn)錄組，所需要的樣品量少，并且可以測(cè)定單細(xì)胞內(nèi)的RNA表達(dá)水平;Pacbio技術(shù)可以直接測(cè)得接近全長(zhǎng)的轉(zhuǎn)錄本，可以有效解決二代技術(shù)拼接較為零碎以及Win10BashonUbuntuonWindows怎么用？（安裝和用法mkdirBiosoftmkdirBiosoft&cd /jmzeng/RNA-seq/RNA-seq-example-GSE81916-two-fastqchisatscdcdmkdirsratoolkit&&cd#長(zhǎng)度：7647376973Mapplication/x-正在保存至“sratoolkit.2.8.2-1-tar-zxvfsratoolkit.2.8.2-1-cdmkdirfastqc&&cd#unzip6.unzipjava，根據(jù)系統(tǒng)提示安裝，一般是sudoaptinstallcdmkdirHISAT&&cd#/pub/infphilo/hisat2/downloads/hisat2-2.1.0-cdhisat2-makeg++gcc，同上用sudoaptinstallrm-f*.hcdmkdirsamtools&&cd#/samtools/samtools/releases/download/1.5/samtools-長(zhǎng)度：4190142(4.0Mapplication/octet-正在保存至samtools-tarxvfjsamtools-1.5.tar.bz2cdsamtools-1.5sudomakesudoapt-getinstallpython-pipinstall--upgradesudoapt-getinstallbuild-essentialpython2.7-devpython-numpypython-#tarxvzfHTSeq-pythonsetup.pyinstall--~/.local/bin/htseq-count--#viexport condainstall-cbiocondacondainstall-cbiocondahtseq=0.7.2condainstall-cbiocondahisat2=2.1.0condainstall-cbioconda10.condainstall-cjfear /p/27670618?utmsource=wechatsession&utmmedium=so AKAP95regulatessplicingthroughscaffoldingRNAsandRNAprocessing#esearchdbsraqueryPRJNA299273|efetchformatruninforuninfo.txt這個(gè)命SRR#catruninfo.txt|cut-f1-d','|grepSRR>~/biosoft/sratoolkit.2.8.2-1-centos_linux64/bin/prefetchoption-filesra.ids數(shù)據(jù)存在/home/shenmy/ncbi/public/sra這個(gè)文件下面，找了半天mkdir/mnt/d/rna_seq/data&&cdperl-lne'$id=substr($_,0,6);print"axelftp://ftp- RR/$id/$_/$_.sra"'SRR_Acc_List.txtbash改成用axel下是ls*.sra|whilereadid;do(/mnt/d/Software/Biosoft/sratoolkit/sratoolkit.2.8.2-ubuntu64/bin/fastq-dump--split-3rmodu-w/ncbinlm /s/LNrQNSHdcdKb25s4mFeqw 了解fastq數(shù)sra文件轉(zhuǎn)換為fastq /mnt/d/Software/Biosoft/sratoolkit/sratoolkit.2.8.2-1-ubuntu64/bin/fastq-dump-Legacy3-filesplittingformate-pairs: biologicalreadssatisfyingdumconditionsareplacedinfiles*_1.fastqand*_2.fastqIfonlyonebiologicalreadispresentitisplacedin*.fastqBiologicalreadsandaboveare也就是說(shuō)如果 Compressoutputusing將這些文件合并成一個(gè)tar文件，然后再使用gzip進(jìn)行壓縮，最后生成的.tar.gz或者.tgz文件就是所謂的“tar壓縮包”或者“tarball”）Compressoutputusing1.1.ls*.sra|whilereadid;do(/mnt/d/Software/Biosoft/sratoolkit/sratoolkit.2.8.2-1-ubuntu64/bin/fastq-dump--split-3$id);done的數(shù)據(jù)是Illumina的雙端4，所以用fastq-dump--split-3命令來(lái)把sra格式數(shù)據(jù)轉(zhuǎn)換為fastq。perl-F'\t'-alneperl-F'\t'-alne"$F[7]\t$F[6]_$F[13]"}'SraRunTable.txt>Rename.txtperl-F'\t'-alne'print"/mnt/d/Software/Biosoft/sratoolkit/sratoolkit.2.8.2-1-ubuntu64/bin/fastq-dump--split-3--gzip-A$F[1]$F[0].sra"'Rename.txt>rmodu-w在此最開(kāi)始之前應(yīng)該做

數(shù)據(jù)傳輸完整性驗(yàn)證1.md5sum*.fastq.gz1.md5sum*.fastq.gz質(zhì)控mkdirmkdir/mnt/d/rna_seq/work&&cdln/mnt/d/rna_seq/data/*mkdirls*.gz|whilereadid;do(/mnt/d/Software/Biosoft/fastqc/FastQC/fastqc/mnt/d/rna_seq/work/1_FastQC_Raw_Data-t-foriin*.zip;dounzip$i;perl|/mnt/d/Software/Biosoft/csvtk/csvtk293cell21.在UCSChg19參考組2.從gencode數(shù)據(jù) 5.在 hg19參 這個(gè)對(duì)新手來(lái)說(shuō)，是一個(gè)很大的坑，hg19、GRCH37、ensembl75這3種 的是同樣的fasta序列， 單位，即NCBI，UCSC及ENSEMBL各自發(fā)布的 樣，比如BGI做的炎黃 組，還有DNA雙螺旋結(jié)構(gòu)提出者 有2016 hg19和hg38，都是UCSC提供的，雖然hg38相比hg19來(lái)說(shuō)，做了很多改進(jìn)，優(yōu)點(diǎn)也不少，但因?yàn)槟壳盀橹购芏嘧⑨屝畔⒍际轻槍?duì)于hg19的坐標(biāo)系統(tǒng)來(lái)的， mkdir/mnt/d/rna_seq/data/reference&&cdmkdir-pgenome/hg19&&cd#nohup nohup tarzvfxcat*.fa>rmGRCh37ist omeReferenceConsortiumHumangenomebuild37.hg19組大小是3G，壓縮后百兆。從gencode數(shù)據(jù)庫(kù)注釋文件，并且用IGV去查看感的的結(jié)gzip-dIGV、tar-zxvfbedtools-cd截圖幾個(gè)的IGV可視化結(jié)grepgrep-w'gene'gencode.v26lift37.annotation.gtf1,4,5grep-w'gene'gencode.v26lift37.annotation.gtf|grep1,4,5>>gene.bedgrep-w'gene'gencode.v26lift37.annotation.gtf1,4,5-w'TP53'|cut--w'KRAS'|cut--w'EGFR'|cut-/mnt/d/Software/Biosoft/bedtools2/bin/bedtoolsigv-igene.bed#perl-alne'{print"goto$F[0]:$F[1]-$F[2]\nsnapshot$F[3].png"}ENSEMBL，NCBI的gtf#axel/pub/grch37/release-3.axelftp://ftp.ensembl. ## 4. 出來(lái)以上里的，還有Bowtie2,BBMap,BWA,CLC,Novoalign,SMALT等，百科hisat2本作業(yè)是比對(duì)到組，所以使用gappedorsplicesmapper，此流程已經(jīng)更新。TopHat首次被已經(jīng)是7年前，STAR的比對(duì)速度是TopHat的50倍，HISAT更是STAR的1.2倍。10HISAT2TopHat2/Bowti2BWT算法，實(shí)現(xiàn)了更快的速度和更少的資源占用，作者TopHat2/Bowti2和HISAT的用戶(hù)轉(zhuǎn)換到HISAT2。11HELP$./hisat2-HISAT2version2.1.0byDaehwanKim hisat2[options]*-x<ht2-idx>{-1<m1>-2<m2>|-U<r>|--sra-acc<SRAaccessionnumber>}[-S<sam>]<ht2-idx>Indexfilenameprefix(minustrailing Fileswith#1mates,pairedwithfilesinCouldbegzip'ed(extension:.gz)orbzip2'ed(extension: Fileswith#2mates,pairedwithfilesinCouldbegzip'ed(extension:.gz)orbzip2'ed(extension: FileswithunpairedCouldbegzip'ed(extension:.gz)orbzip2'ed(extension:<SRAaccession Comma-separatedlistofSRAaccessionnumbers,e.g.-sra-acc FileforSAMoutput(default:stdout)<m1>,<m2>,<r>canbecomma-separatedlists(nowhitespace)andcanspecifiedmanytimes.E.g.'-Ufile1.fq,file2.fq-Ufile3.fq'.20.Options(defaultsin ueryinputfilesareFASTQ.fq/.fastq-- ueryinputfilesareinIllumina'sqseq- queryinputfilesare(multi-)FASTA- queryinputfilesarerawone-sequence-per-- <m1>,<m2>,<r>aresequencesthemselves,not-s/--skip skip <int>reads/pairsintheinput-u/--upto stop <int>reads/pairs(no-5/--trim5 trim<int>basesfrom5'/leftendofreads-3/--trim3 trim<int>basesfrom3'/rightendofreads-- qualitiesarePhred+33-- qualitiesare--int- qualitiesencodedasspace-delimited SRAaccessionID--n-ceil funcformax#non-A/C/G/Tspermittedinaln--ignore- treatallqualityvaluesas30onPhredscale-- donotalignforward(original)versionofread donotalign plementversionofread(off)Spliced--pen-cansplice penaltyforacanonicalsplicesite penaltyforanon-canonicalsplicesite--pen-canintronlen<func> penaltyforlongintrons(G,-8,1)withcanonicalsplicesites--pen-noncanintronlen<func> penaltyforlongintrons(G,-8,1)withnoncanonicalsplicesites--min-intronlen minimumintronlength--max-intronlen umintronlength providealistofknownsplice--novel-splicesite-outfile<path>reportalistofsplice providealistofnovelsplice--no-temp- disabletheuseofsplicesites--no-spliced- disablespliced--rna-strandness specifystrand-specificinformation reportsonlythosealignmentswithinknown reportsalignmentstailoredfortranscript reportsalignmentstailoredspecificallyfor triestoavoidaligningreadstopseudogenes(experimentaloption) disablestemplaengthadjustmentforRNA-seq--mp maxandminpenaltiesformismatch;lowerqual=lower--sp maxandminpenaltiesforsoft-clip;lowerqual=penalty65.--no- nosoft-66.--np penaltyfornon-A/C/G/Tsinread/ref67.--rdg<int>,<int>readgapopen,extendpenalties--rfg<int>,<int>referencegapopen,extendpenalties--score-min<func>minacceptablealignmentscorew/r/tread(L,0.0,--k<int>(default:5)reportupto<int>alnsperreadPaired--I/--minins<int>minimumfragmentlength(0),onlyvalidwith--no-spliced--X/--maxins<int> umfragmentlength(500),onlyvalidwith--no-spliced---fr/--rf/-- -1,-2matesalignfw/rev,rev/fw,fw/fw(----no- suppressunpairedalignmentsforpaired suppressdiscordantalignmentsforpairedreads-t/-- printwall-clocktimetakenbysearch--un writeunpairedreadsthatdidn'talignto--al writeunpairedreadsthatalignedatleastonceto--un-conc writepairsthatdidn'talignconcordantlyto--al-conc writepairsthatalignedconcordantlyatleastonceto(Note:for--un,--al,--un-conc,or--al-conc,add'-gz'totheoptionname,--un-gz<path>,togzipcompressoutput,oradd'-bz2'tobzip2compress--summary- printalignmentsummarytothis printalignmentsummaryinanewstyle,whichismoremachine--- printnothingtostderrexceptserious--met-file<path>sendmetricstofileat<path>--met- sendmetricstostderr--met reportinternalcounters&metricsevery<int>secs--no- supppressheaderlines,i.e.linesstartingwith--no- supppress@SQheader--rg-id setreadgroupid,reflectedin@RGlineandRG:Z:opt--rg add<text>("lab:value")to@RGlineofSAMNote:@RGlineonlyprintedwhen--rg-idis put'*'inSEQandQUALfieldsforsecondary-o/--offrate<int>overrideoffrateofindex;mustbe>=index's-p/--threads<int>numberofalignmentthreadstolaunch-- forceSAMoutputordertomatchorderofinput usememory-mappedI/Oforindex;many'hisat2'scan--qc-filteroutreadsthatarebadaccordingtoQSEQ--seedseedforrandomnumbergenerator --non-deterministicseedrand.gen.arbitrarilyinsteadofusingread--remove-remove'chr'fromreferencenamesin--add-add'chr'toreferencenamesin--printversioninformationand-h/--printthisusagecdcdmkdir-pindex/hisat&&cdwget-cwget-ctarzxvftarxvzf把fastq格式的reads比對(duì)上去得到sam文件，接著用samtools把它轉(zhuǎn)為bam文件，并且接輸出BAM文件12。~/biosoft/HISAT/current/hisat2-p5-x$reference-USRR3589959.fastq-Scontrol_1.sam2>control_1.log~/biosoft/HISAT/current/hisat2-p5-x$reference-USRR3589960.fastq-Scontrol_2.sam2>control_2.log~/biosoft/HISAT/current/hisat2-p5-x$reference-USRR3589961.fastq-SAkap95_1.sam~/biosoft/HISAT/current/hisat2-p5-x$reference-USRR3589962.fastq-SAkap95_2.samls*sam|whilereadid;do(nohupsamtoolssort-n-@5-o${id%%.*}.Nsort.bam$id在上面的基礎(chǔ)進(jìn)行修改，編寫(xiě)bash #!set-uset-set-o11.ls--color=neverHomo*1.fastq.gz|whilereadid;do(~/biosoft/hisat2-2.1.0/hisat2-t-$NUM_THREADS-x$hg19_ref-1$data_path/${id%_*}_1.fastq.gz-$data_path/${id%_*}_2.fastq.gz2>${id%_*}_map.log|samtoolsview--13.ls--color=neverMus*1.fastq.gz|whilereadid;do(~/biosoft/hisat2-2.1.0/hisat2-t-p$NUM_THREADS-x$mm10_ref-1$data_path/${id%_*}_1.fastq.gz-2$data_path/${id%_*}_2.fastq.gz2>${id%_*}_map.log|samtoolsview--15.ls--color=never*.bam|whilereadid;do(samtoolssort--threads$NUM_THREADS$id-16.ls--color=never*_sorted.bam|whilereadid;do(samtoolsindex：bashbamreadsreads比對(duì) 組序列,多重比對(duì)~/biosoft/qualimap_v2.2.1/qualimaprnaseq-Homo_sapiens_Control_293_cell_sorted.bam-outdirHomo_sapiens_Control_293_cell_sorted_QC-pe-s--java-mem-size=60G-gtfcondainstallcondainstallucsc-gtfToGenePred-genePredExt-ignoreGroupsWithoutExons-geneNameAsName2gencode.v26lift37.annotation.gtfgencode.v26lift37.annotation.gpdgzip-djava-jar~/biosoft/picard_2.8.0/picard.jarCollectRnaSeqMetrics bam_stat.py-i載入IGV 看上去很類(lèi)似fastq文件，它也有readread名稱(chēng)mate名稱(chēng)，記錄matepairchromosomeCIGAR 詳情samreadnamereadscountls--color=never*.bam|whilereadid;do(samtoolssort-n--threads40$id-ls--color=never*_sorted_name.bam|whilereadid;do(samtoolsindexreadsname#[E::hts_idx_push]unsorted#samtoolsindex:"Homo_sapiens_AKAP95_KD_miR_12_293_cell_sorted_name.bam"iscorruptedorunsortedsamtoolssort-n--threads30Homo_sapiens_AKAP95_KD_miR_12_293_cell.bam-osamtoolsviewHomo_sapiens_AKAP95_KD_miR_12_293_cell_sorted_name.bam|less-pos排序的文件小些，為什么呢？因?yàn)榘凑瘴恢门判蛳嗨频膬?nèi)容會(huì)排在一起，按照采用默認(rèn)的pos排序。提高文件的壓縮比例。BAMsort縮比提高了，因此排序之后的BAM文件變小了，相對(duì)應(yīng)的SAM文件就是純文本文件，對(duì)SAM文件進(jìn)行排序就不會(huì)改變文件大小。而且由于RNA-seq中由于表達(dá)量的關(guān)系，RNA-seq的數(shù)據(jù)比對(duì)結(jié)果BAM文件使用samtools進(jìn)行sort之后文件壓縮比例變化會(huì)比DNA-seq更甚。另外，samtoolsBAMreads會(huì)被放From:Gainingcomprehensivebiologicalinsightintothetranscriptomebyperformingabroad-spectrumRNA-seq FigureTheRNACocktail ysisprotocol.RNACocktailisacomprehensiveprotocolofRNA-seq ysis.ThefiguresummarizesthewidelyusedapproachesforthekeystepsoverthebroadspectrumofRNA-seq ysisandalsosuccinctlycapturesthepossibleworkflowsonecanuseto yseRNA-seqdatabedtoolsusageusage:htseq-count[options]alignment_file<alignment_file>：比對(duì)到組后得到的SAM文件（SAMtools包含一些perl可以將大多數(shù)的比對(duì)文件轉(zhuǎn)換成SAM格式），注意組map時(shí)一定要用支持剪接的比對(duì)（splicing-awarealigner）進(jìn)行比對(duì)如TopHat.HTSeq-count需要用到SAM格式中的CIGAR區(qū)域的信息。<gff_file>:包含單位信息的gff/GTF文件（gff文件格式），大多數(shù)情況下就是指注釋文件 由GTF文件其實(shí)就是gff文件格式的變形，在這里同樣可以傳入GTFThisscripttakesoneormorealignmentfilesinSAM/BAMformatandafeaturefileinGFFformatandcalculatesforeachfeaturethenumberofreadsmap toit.Seeforpositional PathtotheSAM/BAMfilescontainingthemapped If'iselectd,adfrmstndardinut想通過(guò)準(zhǔn)輸入來(lái)傳入 組map 到SAM文，用–替換<aligmetfile>即可 Pathtothefilecontainingthefeaturesoptional-h,-- showthishelpmessageand-f{sam,bam},--formattypeof<alignment_file>data,either'sam'or(default:指定輸入文件的格式，<format可以是samfortextSAMfiles)bam(forbinaryBAMfiles)-r{pos,name},--order'pos'or'name'.Sortingorderof(default:name).**Paired-endsequencingdatamustbesortedeitherbypositionorbyreadname,andthesortingordermustbespecified.**Ignoredforsingle-enddata.如果你是雙端，必須要對(duì)SAM進(jìn)行排序（單端可不必排序，但這里我也對(duì)單端結(jié)果排序已減少內(nèi)存消耗并提高效率），對(duì)readname或 排序皆可（這里我按readname排序，因?yàn)橥ㄟ^(guò)位置排序我遇到過(guò)錯(cuò)誤）。具體需要通過(guò)-r參數(shù)指定，所以排序請(qǐng)?jiān)斠?jiàn)參數(shù)-r，在sort時(shí)用-n則不用修改，默認(rèn)的排序則修改成pos。When<alignment_file>ispairedendsortedbyposition,allowonlysomanyreadstostayinmemoryuntilthematesarefound(raisingthisnumberwillusemorememory).Hasnoeffectforsingleendorpairedendsortedbyname-s{yes,no,reverse},--strandedwhetherthedataisfromastrand-specificassay.Specify'yes','no',or'reverse'(default:yes).'reverse'means'yes'withreversedstrand是否鏈特異性，如果不是修改成-aMINAQUAL,--minaqualskipallreadswithalignmentqualitylowerthanthegivenminimumvalue(default:指定一個(gè)最低read -tFEATURETYPE,--typefeaturetype(3rdcolumninGFFfile)tobeused,allfeaturesofothertypeareignored(default,suitableforEnsemblGTFfiles:exon)指定最小計(jì)數(shù)單位類(lèi)型（gff文件的第3列中的類(lèi)型如：exon）,指定后其他單位類(lèi)型將被忽略（默認(rèn)情況下，對(duì)于rna-seq分析采用EnsemblGTF： -iIDATTR,--idattrGFFattributetobeusedasfeatureID(default,suitableforEnsemblGTFfiles:最終的計(jì)數(shù)單位，一般為。默認(rèn)為：gene_id需要改成gene_name，可以直接識(shí)別，如果是gene_id--additional-attrADDITIONAL_ATTR[ADDITIONAL_ATTRAdditionalfeatureattributes(default:none,suitableforEnsemblGTFfiles:-m{union,intersection-strict,intersection-nonempty},--mode{union,intersection-modetohandlereads morethanone(choices:union,intersection-strict,intersection-nonempty;default:判斷一個(gè)reads屬于某個(gè)的模型，用來(lái)判斷統(tǒng)計(jì)reads的時(shí)候?qū)σ恍┍容^特殊的reads定義是否計(jì)入<mode>unionintersection-strict、intersectionnonempty（--nonuniqueWhethertoscorereadsthatarenotuniquelyalignedorambiguouslyassignedto-oSAMOUTS[SAMOUTS...],--samoutSAMOUTS[SAMOUTSwriteoutallSAMalignmentrecordsintoanoutputSAMfilecalledSAMOUT,annotatingeachlinewithitsfeatureassignment(asanoptionalfieldwithtag'XF')出所有algnment的reads叫<samout>的sam件中通過(guò)個(gè)可選的sam ‘XF來(lái)標(biāo)注一對(duì)應(yīng)的位計(jì)數(shù)，以設(shè)-q, s