在將原始數(shù)據(jù)和參考基因組數(shù)據(jù)處理好以后，就開始開始比對分析了
比對所用到的參考基因組的索引文件和基因組注釋文件都存放在hg19文件夾
將sra文件解壓至Fastq文件夾
主要步驟有下列幾步
1用Ttimmomatic對fastq數(shù)據(jù)去接頭

mkdir trim_out #創(chuàng)建一個存放用Trimmomatic去接頭的輸出文件夾
for i in Fastq/*fastq
do
        echo $i
        a=$(echo $i | cut -d "/" -f5 | cut -d "_" -f1)
        java -jar ~/Biosofts/Trimmomatic-0.38/trimmomatic-0.38.jar SE -phred33 $i ~/trim_out/$a-clean.fastq ILLUMINACLIP:/home/qiujunhui/Biosofts/Trimmomatic-0.38/adapters/TruSeq2-SE.fa:2:30:10
done

2.序列比對用到tophat2軟件，使用tophat軟件的優(yōu)點(diǎn)在于tophat2在將待測序列與參考基因組比對后，會直接生成bam文件，生成的bam文件直接可以給cufflinks構(gòu)建轉(zhuǎn)錄本

mkdir tophat_out #創(chuàng)建一個用tophat比對的輸出文件夾
for x in ~/trim_out/*clean.fastq
do
          echo $x
          b=$(echo $x | cut -d "/" -f5 | cut -d "-" -f1)
          mkdir ~/tophat_out/$b
          tophat -p 4 -o ~/tophat_out/$b ~/hg19/genome $x
done            

-p 指定線程
-o 指定輸出目錄
hg19/genome為bowtie2官網(wǎng)上人類基因組的索引文件

2.利用cufflinks構(gòu)建轉(zhuǎn)錄本

#mkdir cufflinks_out 創(chuàng)建一個用cufflinks構(gòu)建轉(zhuǎn)錄本的輸出文件夾
cd cufflinks_out
mkdir SRR957677
mkdir SRR957678
mkdir SRR957679
mkdir SRR957680
cd ..
cufflinks -g ~/hg19/gencode.v29lift37.annotation.gtf -o ~/cufflinks_out/SRR957677 -p 4 ~/tophat_out/SRR957677/accepted_hits.bam
cufflinks -g ~/hg19/gencode.v29lift37.annotation.gtf -o ~/cufflinks_out/SRR957678 -p 4 ~/tophat_out/SRR957678/accepted_hits.bam
 cufflinks -g ~/hg19/gencode.v29lift37.annotation.gtf -o ~/cufflinks_out/SRR957679 -p 4 ~/tophat_out/SRR957679/accepted_hits.bam
cufflinks -g ~/hg19/gencode.v29lift37.annotation.gtf -o ~/cufflinks_out/SRR957680 -p 4 ~/tophat_out/SRR957680/accepted_hits.bam

3.利用cuffmerge合并轉(zhuǎn)錄組
為了比較不同樣品的差異，需要將實(shí)驗(yàn)組和對照組的轉(zhuǎn)錄組合并起來，cuffmerge不僅能用來合并，兩個或多個轉(zhuǎn)錄組，還能把注釋過后的基因組信息也合并起來，從而找到新的基因可變剪切提高合并轉(zhuǎn)錄組的質(zhì)量

mkdir cuffmerge_out #創(chuàng)建一個存放用cuffmerge合并轉(zhuǎn)錄本的輸出文件夾
#把上一操作得到的transcripts.gtf的絕對路徑寫到一個文本文件里面
vi accepted.txt
~/home/qiujunhui/cufflinks_out/SRR957677/transcripts.gtf
~/home/qiujunhui/cufflinks_out/SRR957678/transcripts.gtf
~/home/qiujunhui/cufflinks_out/SRR957679/transcripts.gtf
~/home/qiujunhui/cufflinks_out/SRR957680/transcripts.gtf
#保存
cuffmerge -g ~/hg19/gencode.v29lift37.annotation.gtf -o ~/cuffmerge_out -p 4 accpeted.txt

4.利用cuffdiff進(jìn)行基因表達(dá)差異分析

mkdir cuffdiff _out #創(chuàng)建一個存放用cuffdiff進(jìn)行基因表達(dá)差異分析的輸出文件夾
cuffdiff -o ~/cuffdiff_out -L lable1,lable2,lable3,lable4 -p 4 -u ~/cuffmerge_out/merged.gtf ~/tophat_out/SRR957677/accepted_hits.bam  ~/tophat_out/SRR957677/accepted_hits.bam 
~/tophat_out/SRR957677/accepted_hits.bam ~/tophat_out/SRR957677/accepted_hits.bam

merged.gtf為上一步用cuffmerge合并的轉(zhuǎn)錄組
-L 后面是bam文件的標(biāo)簽，有幾個bam文件就取幾個標(biāo)簽，我只有四個bam文件所以只有四個標(biāo)簽

cuffdiff輸出文件比較多，它會對每個基因，每個轉(zhuǎn)錄片段，每個編碼序列，每個基因的不同剪切體進(jìn)行FPKM，個數(shù)和樣本間差異進(jìn)行分析，最后生成機(jī)組不同的文件，按照不同的需求，就可以往下分析了
cuffdiff計(jì)算每個樣本中的轉(zhuǎn)錄本，初始轉(zhuǎn)錄本和基因FPKM
1.traisoforms.fpkm_tracking 轉(zhuǎn)錄組的FPKM

gens.fpkm_tracking 基因的fpkm

cds fpkm_tracking 編碼序列的fpkm

tss_groups.fpkm_tracking 原始轉(zhuǎn)錄組的FPKM

2.Count tracking files

評估每個樣本中來自每個 transcript, primary transcript,和 gene的fragment數(shù)目

3.Read group tracking

計(jì)算在每個repulate中每個transcript， primary transcript和gene的表達(dá)量和frage數(shù)目

4.Differential expression test

對于splicing transcript， primary transcripts, genes,

and coding sequences.樣本之間的表達(dá)差異檢驗(yàn)。