前一段時間跟著孟浩巍大神的視頻學(xué)習(xí)动雹，在自己的小破筆記本上還是跑完了整個RNAseq差異表達(dá)的分析流程（ tophat2 + cufflink + cuffdiff ）雖然這個流程比較老了鳖枕，現(xiàn)在做分析一般使用的都是 HTseq + DESeq2 等其他的流程祝懂，但是作為入門的知識還是比較容易理解，這篇文章先更一下流程涩拙，后面會再更一篇 安裝 Linux 子系統(tǒng)（已更新）飘蚯，安裝 anconda 和一些分析軟件（已更新）的流程，湊一個真正完整的入門生信的操作流程艳吠。
火山圖麦备、熱圖在 R 中可視化部分已更新

我的電腦配置，真的是戰(zhàn)五渣昭娩。

電腦配置

但還是在一天內(nèi)跑完了整個流程

運行環(huán)境python2.7

原始數(shù)據(jù)如下：

原始文件

包括四個文件：

• bowtie2_hg19 index 文件（這里已經(jīng)提前使用bowtie2建立好了index文件可以直接使用）
• raw_data illumina 雙端測序原始文件（對照組和實驗組各兩個凛篙，就是八條測序文件）
• rRNA rRNA index 序列文件（用于去除 rRNA 的影響）
• RefSeq_genes_hg19.gtf 基因組注釋文件
  鏈接：https://pan.baidu.com/s/1Q4SwosKaG16-aYA74SaASQ
  提取碼：gss3

分析流程

• RNA-seq的原始數(shù)據(jù)（raw data）的質(zhì)量評估
• raw data的過濾和清除不可信數(shù)據(jù)（clean reads）
• reads回帖基因組和轉(zhuǎn)錄組（alignment）
• 計數(shù)（count ）
• 基因差異分析（Gene DE）
• 數(shù)據(jù)的下游分析（這次先不做這個，下次會單獨寫）

下面開始正式分析

1栏渺、fastqc質(zhì)控

mkdir fastqc_result.raw     #(建立輸出文件夾)                                                                       
fastqc -q -t 3 -o ../fastqc_result.raw/ *.fq.gz &    #（使用fastqc質(zhì)控軟件呛梆，對所有raw_data進(jìn)行質(zhì)控檢測）

2、multiqc整合質(zhì)控文件（因為得到很多的質(zhì)控檢測結(jié)果磕诊，需要整合）

multiqc .  #(這一步就是對 fastqc_reault 文件夾下所有文件進(jìn)行整合)

整合后文件

3填物、根據(jù)質(zhì)控結(jié)果纹腌，使用fastx_tolltik去除不良序列

zcat hela_ctrl_rep1_R1.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_rep1_R1.fq.gz &             
zcat hela_ctrl_rep2_R1.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_rep2_R1.fq.gz &             
zcat hela_ctrl_rep2_R2.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_rep2_R2.fq.gz &             
zcat hela_ctrl_rep1_R2.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_rep1_R2.fq.gz &             
zcat hela_treat_rep1_R1.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_t_rep1_R1.fq.gz &          
zcat hela_treat_rep1_R2.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_t_rep1_R2.fq.gz &          
zcat hela_treat_rep2_R1.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_t_rep2_R1.fq.gz &          
zcat hela_treat_rep2_R2.fq.gz  | fastx_trimmer -f 11 -l 140 -z -o out_t_rep2_R2.fq.gz &          

因為當(dāng)時我還不會寫 shell 腳本，正則表達(dá)式也不懂滞磺，就一個一個寫了壶笼，但是 shell 才是生產(chǎn)力啊啊啊啊,這一步也可以放后臺的，當(dāng)時為了看結(jié)果就一個一個跑了

zcat解壓縮雁刷，文件名，fastx_trimmer -f x -l y 保留從x-y的序列 -z壓縮命令 -o輸出結(jié)果 &后臺運行
trimmer保礼，可以使所有序列長度一致

4沛励、cutadaptor去掉read兩端的adaptor

nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_c_rep1_R1.fq.gz -p cut_out_c_rep1_R2.fq.gz out_c_rep1_R1.fq.gz out_c_rep1_R2.fq.gz &       
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_c_rep2_R1.fq.gz -p cut_out_c_rep2_R2.fq.gz out_c_rep2_R1.fq.gz out_c_rep2_R2.fq.gz &          
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_t_rep1_R1.fq.gz -p cut_out_t_rep1_R2.fq.gz out_t_rep1_R1.fq.gz out_t_rep1_R2.fq.gz &          
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_t_rep2_R1.fq.gz -p cut_out_t_rep2_R2.fq.gz out_t_rep2_R1.fq.gz out_t_rep2_R2.fq.gz &          

--times 1一條序列只去一次Adaptor；
-e 0.1在匹配時可以有10%的錯誤率炮障；
-o 3 adaptor序列必須和測序序列有3個堿基以上的overlap才可以目派；
--quality-cutoff常用6；
-m 50如果處理之后低于50的話就扔掉序列胁赢，短序列測序質(zhì)量可能不是很好企蹭；
-a和-A是 Illumina 常用的通用引物，之所以輸入兩個智末，是因為是一個雙端測序的結(jié)果谅摄，需要對兩個文件內(nèi)容進(jìn)行分別去除，-a對應(yīng)Reads1系馆，-A對應(yīng)reads2
-o 上一步輸出fastq1 -p 上一步輸出fastq2
> 最后是寫入log文件
//其中nohup：不掛斷地運行命令送漠。
& 就是放后臺
//2>1 是傳遞給shell腳本的第一個參數(shù)；

5由蘑、利用bowtie2將reads比對到rRNA上闽寡，除去rRNA的影響

nohup bowtie2 -x $rRNA_index -1 cut_out_c_rep1_R1.fq.gz  -2 cut_out_c_rep1_R2.fq.gz  -S sam_out_c_rep1_bowtie -p 2 --un-conc-gz fastq_unmap_c_rep1_R1 &
nohup bowtie2 -x $rRNA_index -1 cut_out_c_rep2_R1.fq.gz  -2 cut_out_c_rep2_R2.fq.gz  -S sam_out_c_rep2_bowtie -p 2 --un-conc-gz fastq_unmap_c_rep2_R1 &                                                  
nohup bowtie2 -x $rRNA_index -1 cut_out_t_rep2_R1.fq.gz  -2 cut_out_t_rep2_R2.fq.gz  -S sam_out_t_rep2_bowtie -p 2 --un-conc-gz fastq_unmap_t_rep2_R1 &                                                  
nohup bowtie2 -x $rRNA_index -1 cut_out_t_rep1_R1.fq.gz  -2 cut_out_t_rep1_R2.fq.gz  -S sam_out_t_rep1_bowtie -p 2 --un-conc-gz fastq_unmap_t_rep1_R1 &  

（這里要先為rRNA進(jìn)行index建庫，如果有別人建好的庫可以直接下載使用）

 rRNA_index=/mnt/c/Users/wt/Desktop/test_data/rnaseq_test_date/rRNA/bowtie2_rRNA_human

一般通過map到rRNA中的比例來衡量建庫的質(zhì)量尼酿。一般的要求rRNA的比例不超過10%。
-x 對應(yīng)rRNA的索引序列裳擎；
-1涎永，-2 是剛輸出的reads1和reads2土辩；
-S 是比對結(jié)果的輸出文件启涯；
-p 2 使用2個核心去運算（p就是processor吧！）蒸殿；
--un-conc-gz 輸出比對不上的結(jié)果；(比對不上的才是需要的)

輸出結(jié)果如下

輸出結(jié)果

其實還應(yīng)當(dāng)將log文件輸出霞玄，用于查看運行過程喊暖，如果報錯容易查看進(jìn)行處理咨跌，但是我第一次做的時候輸出的都是空白，也不知道哪里有問題，遍膜。這里sam_out文件有點問題弛说，雖然沒有用，本應(yīng)當(dāng)輸出的是比對上的比例翰意，是一個log文件木人，后面會再看看這里怎么回事。

到這里才算質(zhì)控做完冀偶！

6虎囚、使用 tophat2 將過濾掉 rRNA 的 reads 比對到 ref（參考）基因組上

（如果mRNA直接比對到人的DNA上，可能會出現(xiàn)問題蔫磨，有可能跨越了一個內(nèi)含子，tophat2考慮了這個問題圃伶，它將reads根據(jù)注釋文件分開成短序列堤如，重新比對；）
需要先建庫（這里別人建好了）

 hg19_index=/mnt/c/Users/wt/Desktop/test_data/rnaseq_test_date/bowtie2_hg19/hg19_only_chromosome  

nohup tophat2 -p 2 -o top_out1 $hg19_index fastq_unmap_c_rep1_R1.fq.1.1.gz fastq_unmap_c_rep1_R1.fq.1.2.gz &
nohup tophat2 -p 2 -o top_out2 $hg19_index fastq_unmap_c_rep2_R1.fq.2.1.gz fastq_unmap_c_rep2_R1.fq.2.2.gz & 
nohup tophat2 -p 1 -o top_out3 $hg19_index fastq_unmap_t_rep1_R1.fq.1.1.gz fastq_unmap_t_rep1_R1.fq.1.2.gz & 
nohup tophat2 -p 1 -o top_out4 $hg19_index fastq_unmap_t_rep2_R1.fq.2.1.gz fastq_unmap_t_rep2_R1.fq.2.2.gz & 

-o top_out1窒朋， 結(jié)果輸出到這個文件夾中
hg19 是人的第19個基因組版本搀罢，常用的包括hg19和hg38，hg38會取代hg38侥猩，hg19包含的信息比較豐富
所使用序列是上一步未比對上的序列unmap(也就是我們所需要的質(zhì)控后的結(jié)果)

輸出結(jié)果如下

tophat2輸出結(jié)果

.bam 是最終的比對結(jié)果;
.txt 是比對中的總結(jié)情況榔至；
3個bed文件，軟件檢測出來的 deletions 缺失欺劳， insertions 插入唧取， junctions 區(qū)域
.info 有一些reads沒有直接 map 到連續(xù)的 DNA 基因組上，需要切一些reads,加工 reads 的過程文件保存在 info 里划提；
unmapped.bam 是沒有 map 上去枫弟，一層層都沒有比對到的，可能是基因組上未注釋過的鹏往、測序問題等淡诗。

7、使用cuffliks對表達(dá)量進(jìn)行評估(這一步在正常情況下沒什么用)

cufflinks -o cufflink_out1 -p 4 -G ../RefSeq_genes_hg19.gtf top_out1/accepted_hits.bam                                                                                                   
cufflinks -o cufflink_out2 -p 4 -G ../RefSeq_genes_hg19.gtf top_out2/accepted_hits.bam     
cufflinks -o cufflink_out3 -p 4 -G ../RefSeq_genes_hg19.gtf top_out3/accepted_hits.bam      
cufflinks -o cufflink_out4 -p 4 -G ../RefSeq_genes_hg19.gtf top_out4/accepted_hits.bam   

-G 轉(zhuǎn)錄組的參考文件

cufflinks輸出結(jié)果如下：

cufflinks 輸出結(jié)果

gene.fpkm gene的 fpkm 計算結(jié)果（基因表達(dá)量）
isoforms.fpkm isoforms (可以理解為 gene 的各個外顯子)的 fpkm 計算結(jié)果（轉(zhuǎn)錄本表達(dá)量）
skipped.gtf 跳過的基因的轉(zhuǎn)錄本信息
transcripts.gtf 轉(zhuǎn)錄組的gtf,該文件包含Cufflinks的組裝結(jié)果isoforms
fpkm 是衡量基因表達(dá)量的數(shù)值伊履，一個基因有不同的內(nèi)含子和外顯子韩容，不同的外顯子之間可以形成不同的轉(zhuǎn)錄本，每一個轉(zhuǎn)錄本可以翻譯成不同的蛋白唐瀑，這些蛋白互相之間就是 isoforms（亞型）群凶，對于不同的轉(zhuǎn)錄本來說基因有一個表達(dá)量，這就是基因的 fpkm 和 isoform 的 fpkm 哄辣。

8座掘、cuffdiff計算基因表達(dá)差異

ctrl_bam=top_out1/accepted_hits.bam,top_out2/accepted_hits.bam
treat_bam=top_out3/accepted_hits.bam,top_out4/accepted_hits.bam
label=hela_ctrl,hela_treat
cuffdiff -o diff_out1 -p 7 --labels $label --min-reps-for-js-test 2 ../RefSeq_genes_hg19.gtf $ctrl_bam $treat_bam

--lables 是文件的輸入次序递惋，如上 label=hela.ctrl,hela_treat；
--min 每個 treat 里有幾個 repeat 溢陪，上邊 ctrl_bam 是兩個萍虽，要和 treat_bam 數(shù)量一致且>=2（最小重復(fù)）
然后比對到參考轉(zhuǎn)錄本上

運行結(jié)果如下：

cuffdiff 結(jié)果