簡(jiǎn)單地記錄一下我的一次RNA-seq實(shí)戰(zhàn)，由于我們實(shí)驗(yàn)室的服務(wù)器沒(méi)有聯(lián)網(wǎng)，我之前一直在自己的PC上，安裝虛擬機(jī)或者用OS系統(tǒng)跑數(shù)據(jù)，出現(xiàn)許多問(wèn)題，而且速度很慢，后來(lái)詢問(wèn)了洲更大神，直接按照biostar handbook 來(lái)設(shè)置我的電腦Mac，但是我又遇到了很多問(wèn)題，不過(guò)biostar handbook 確實(shí)很不錯(cuò)，我后續(xù)的學(xué)習(xí)應(yīng)該要按照這本書(shū)來(lái)，但是這次時(shí)間比較急，要趕緊先跑一遍，學(xué)會(huì)大致流程再說(shuō)。那我就跌跌撞撞根據(jù)百度，文獻(xiàn)，問(wèn)大神的模式開(kāi)始學(xué)習(xí)了。
整個(gè)學(xué)習(xí)過(guò)程是完全沒(méi)有邏輯和順序的，但是我現(xiàn)在是整理好的順序，這樣后面的人來(lái)看操作起來(lái)就很方便了。
1安裝軟件的過(guò)程
先要安裝這些軟件，這些軟件直接在百度搜索，然后去官網(wǎng)下載軟件包，Linux 下直接解壓可以用的就可以，建議不要源碼編譯，我源碼編譯了一下，但是服務(wù)器沒(méi)有聯(lián)網(wǎng)，缺少lib 庫(kù) 反正沒(méi)有直接下載軟件包來(lái)的方便。
fastqc?
fastx_toolkit
trimmomatic?
bowtie2?
tophat?
cufflinks
1傳輸文件ssh
scp -P 22 -r 要傳輸?shù)奈募窂?will@10.10.10.10 你要保存的文件路徑
2解壓.gz .tar.gz
解壓這些文件，解壓不同類型的方法不一樣，具體可以參考這篇文章：
http://www.cnblogs.com/qq78292959/archive/2011/07/06/2099427.html
*.tar 用 tar -xvf 解壓
*.gz 用 gzip -d或者gunzip 解壓
*.tar.gz和*.tgz 用 tar -xzf 解壓
*.bz2 用 bzip2 -d或者用bunzip2 解壓
3 把文件PATH修改
? 1 cd 回到主目錄
? 2 ls -all 查看所有的文件 找到 ?.bash_profile ( 可能不同系統(tǒng)命名不一樣）
? 3 ?vim .bash_profile
增加 exportPATH=$PATH:/public/home/你的用戶名/rna-seq/FastQC/
4 FASTqc使用指南檢測(cè)測(cè)序的質(zhì)量
這個(gè)看我之前寫(xiě)的一篇文章專門(mén)記錄fastqc 以及遇到的問(wèn)題 或者參考這篇https://zhuanlan.zhihu.com/p/20731723
主要命令如下：
?1 主文件下 chmod 755 fastqc
? 2 fastqc -o result -t 8 gene_1.fq
然后會(huì)生成一個(gè)報(bào)告文件，具體分析可以看之前的文章
5 Trimmomatic處理reads 可以參考
http://www.itdecent.cn/p/36891a89ed6e
http://blog.sina.com.cn/s/blog_4b91a9e50101ock4.html
主要就是java 的程序
java -jar trimmomatic-0.36.jar PE? -threads20 -phred33 -trimlog logfile?reads1.fq reads2.fq ?reads1.clean.fq reads1.unpaired.fq reads2.clean.fq reads2.unpaired.fq ILLUMINACLIP:/path/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10\ LEADING:3TRAILING:3SLIDINGWINDOW:4:15MINLEN:50
PE/SE 設(shè)定對(duì)Paired-End或Single-End的reads進(jìn)行處理，其輸入和輸出參數(shù)稍有不一樣。 -threads 設(shè)置多線程運(yùn)行數(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 Perform a sliding window trimming。Windows的size是4個(gè)堿基，其平均堿基 質(zhì)量小于15，則切除。 MINLEN:50 最小的reads長(zhǎng)度 CROP: 保留reads到指定的長(zhǎng)度 HEADCROP: 在reads的首端切除指定的長(zhǎng)度 TOPHRED33 將堿基質(zhì)量轉(zhuǎn)換為pred33格式 TOPHRED64 將堿基質(zhì)量轉(zhuǎn)換為pred64格式
會(huì)出現(xiàn)4個(gè)文件
http://blog.sciencenet.cn/blog-2458445-933838.html
paired 的一對(duì) unpaired 的一對(duì)
用biotie2 處理paired 的一對(duì)
其實(shí)我還是有點(diǎn)不理解，下次理解了再補(bǔ)上這里的
6 bowtie2 （map的常用工具 ）（短序列拼接至模版基因組）
好了現(xiàn)在reads處理好了，接下來(lái)要map了
?1 先要建立參考基因組的索引 對(duì)參考基因組fasta 格式文件建立index ( 去網(wǎng)上下載參考基因組的fa文件）
bowtie2-build genome.fa genome?
然后你會(huì)看到 genome.1.bt2 4個(gè) 還有 genome.rev.1.bt2 2 個(gè)
7 tophat通過(guò)調(diào)用bowtie2來(lái)完成map
tophat -p 30 -G? annotation.gff3 -o qc-thout genome qc-1.fq qc-2.fq
( annotation.gff3 為基因的注釋信息 genome為基因組的序列，qc-1.fq qc-2.fq為質(zhì)控后的雙端測(cè)序數(shù)據(jù)）
出現(xiàn)錯(cuò)誤一 （注意啦 注意啦）
Error: Could not find Bowtie 2 index files (genome.*.bt2)
這個(gè)是網(wǎng)上找的，和我基本一樣，我也是按照網(wǎng)上攻略輸入的語(yǔ)句，但是就是出現(xiàn)問(wèn)題。然后我搗鼓了好久，去了https://www.biostars.org/p/203950/這個(gè)biostar 網(wǎng)站，還是不錯(cuò)的，不過(guò)里面的英文有些理解不了，要參考一下中文的. 解決如下：
1 最好把要跑的文件和參考基因組的文件都放在一個(gè)文件夾下面
2 (我的問(wèn)題主要是這個(gè)） basename 很重要， 一定要把basename 改成和基因的注釋文件一樣的開(kāi)頭，不能下載下來(lái)是什么就用什么
3 我把基因的注釋文件后面的 .fa拿掉了，發(fā)現(xiàn)才可以跑起來(lái)，不知道為什么
出現(xiàn)錯(cuò)誤二 （ 注意啦注意啦）
跑了大約1個(gè)小時(shí)的時(shí)候突然斷了
Tophat Error "segment-based junction search failed"
這是咋搞得呢，我查了查https://www.biostars.org/p/146557/
這篇英文說(shuō)的還是挺清楚的 我線程用的太多了 我用了 80.....
好吧我把 -P 調(diào)節(jié)到30 就解決了
OK 跑完出現(xiàn)了幾個(gè)文件：
accepted_hits.bam? align_summary.txt? deletions.bed? insertions.bed? junctions.bed? logs? prep_reads.info? unmapped.bam
參考一下這個(gè)解析很清楚：http://www.360doc.com/content/17/0802/19/45954458_676165658.shtml
8 cufflinks查看樣本間基因的差異表達(dá)計(jì)算
主要用于基因表達(dá)量的計(jì)算和差異表達(dá)基因的尋找。
Cufflinks程序主要根據(jù)Tophat的比對(duì)結(jié)果，依托或不依托于參考基因組的GTF注釋文件，計(jì)算出(各個(gè)gene的)isoform的FPKM值，并給出trascripts.gtf注釋結(jié)果(組裝出轉(zhuǎn)錄組)。
1 ?用 cufflinks 對(duì)Tophat比對(duì)的每個(gè)結(jié)果（bam文件）進(jìn)行組裝
? ? ?cufflinks -p 30 -u -g annotation.gff3 --library-type fr-unstranded - o S3 accepted_hits.bam?
出現(xiàn)下面幾個(gè)文件
genes.fpkm_tracking? isoforms.fpkm_tracking? skipped.gtf? transcripts.gtf
2 將 Cufflinks 組裝出來(lái)的轉(zhuǎn)錄本注釋文件（gtf）結(jié)尾的 路徑匯總到文件 assemblies.txt
然后用cuffmerge將各個(gè)樣本的組裝結(jié)果與下載的基因注釋文件整合在一起，創(chuàng)建一個(gè)整體的注釋信息：
cuffmerge -p 30 -o merge_out -g annotation.gff3 -s genome.fa assemblies.txt
注意這個(gè)assemblies.txt文件要自己建立，然后寫(xiě)入路徑
3 最后啦啦啦
用cuffdiff -o diff-out -b genome.fa -p 40 -L S1,G1 -u merged.gtf accepted_hits.bam
9未完待續(xù)后面用R語(yǔ)言來(lái)做可視化
后記：
這個(gè)流程現(xiàn)在看看是很簡(jiǎn)單，但是從一點(diǎn)也不會(huì)到跑下來(lái)遇到的問(wèn)題都是很困難的，特別是沒(méi)有人帶你，還是要感謝各路大神，各種論壇。 說(shuō)不定這個(gè)流程在我服務(wù)器上跑是可以的，但是別人的就不行，所以多多交流哈
