使用的数据,双端测序数据clean data 1 & clean data 2
1、bowtie建库
bowtie2-build genome.fa genome_index (*.fa代表 要建库的文件名字,genome_index代表建库后的名字)
2、将read比对到基因组中
bowtie -q -S -t -p 20 -m 1 --best --strata genome_index input.fastq out.sam
-q 指输入的测序数据是Fastq格式
-S 输出的结果文件格式
-p 使用的线程数
-m 只保留有一个alignment的read
3、查看比对结果
samtools view -S -b sample.sam | samtools flagstat -
但为了使用后面的cufflink套件,一般使用tophat做序列比对
2、使用Tophat做序列比对
依赖软件有bowtie、samtools,不同版本的Tophat对bowtie的要求不同
tophat -p 8 -G annotation.gtf -o tophat_output bowtie_index pair_end-1 pair_end-2
-p 使用的线程数
-G 使用的注释文件
-o 输出文件夹
注意 1、bowtie_index的名字尽量和基因组的名字(不包括后面的.fa)一样
2、此处使用的基因组和注释文件必去相互对应,即基因组各个染色体的名字必须与注释文件的第一列对应
<ATTENTION>
Mapping的最后一步是去除map到基因组中多于一处的序列,如果出现好几个序列都map在完全相同的一个区段,那么就应该只保留一个这样的序列,所以,只保留匹配最高的那一个。而且这样的序列占很大一部分,这步也很简单,samtools里的rmdup可以轻松解决:
samtools rmdup -s input.bam output.bam
3、使用cufflink套件计算表达量和差异表达基因
Cufflinks是一套拼接转录本,计算表达量,计算差异表达的工具。尽可能拼接处最优可能的转录本的结构,并且估计它的表达量。
cufflinks -p 10 -u -g/G annotation.gtf -o Cufflink1_R1_output tophat_output/accepted_hits.bam
-p 使用的线程数
-u 告诉Cufflinks用更准确的方法去处理贴到多个位点上的reads,如果没有-u,Cufflinks只会把这些reads简单地平均分配。比如一个read贴到了10个位置,那么每个位置分得十分之一。加-u后会先进行平均分配,然后按照这10个位置各自的表达量,计算read被分配到每个位置的概率。
-G是提供一个注释文件,并且告诉Cufflinks不要去拼接新的转录本,只能用注释文件里提供的转录本。
-g 也是提供一个注释文件,但是Cufflinks会在这些已知转录本的指导下,拼接新的转录本
-o 输出的文件夹
最后加的bam文件是Tophat运行结束产生的结果
4、cuffmerge整合所有组装好的转录本
当我们使用Cufflinks处理多个数据之后,我们需要将其转录本数据整合为一个全面的转录本集合,Cuffmerge是一个将Cufflinks生成的gtf文件融合为一个更加全面的转录本注释结果的工具。
cuffmerge -g annotation.gtf -s genome.fa -p 10 -o cuffmerge_output assembles.txt
-g 参考注释文件
-s 基因组文件
-o 输出文件夹
-p 使用的线程数
assembles.txt 提供组装的转录本位置
cat > assembles.txt
~/cufflink/C1_R1_output/transcripts.gtf
~/cufflink/C1_R2_output/transcripts.gtf
~/cufflink/C1_R3_output/transcripts.gtf
5、cuffdiff计算差异表达的基因
计算的简单原理在于测序深度和外显子长度一定时,Read的数量与对应的转录本数量成正比。通过对Reads进行计数计算转录本的表达量。同时cuffdiff可以计算不同条件下转录本表达水平的显著性差异。
cuffdiff -o cuffdiff_output -p 10 -L Control_1,Control_2 -u cuffmerge_output/merged.gtf \ ./tophat_output_C1_R1/accepted_hits.bam,./tophat_output_C1_R2/accepted_hits.bam \ ./tophat_output_C2_R1/accepted_hits.bam,./tophat_output_C2_R2/accepted_hits.bam
-p 使用的线程数
-u 同cuffmerge使用的-u参数
-L 为每个样品标上名称
接下来Cuffmerge产生的gtf文件,Cuffdiff需要它提供的注释进行初始转录产物和可变剪切等定量分析。最后是TopHat产生的bam文件,如果一个样品中有多个实验重复,那么我们需要提供哦呢bam文件列表,文件名之间以逗号隔开。
运行之后,cuffdiff输出的文件在diff_out目录之下。其中包括一些按类别统计的表达水平结果,如果有相同的转录起始位点,或具有相同的编码区的转录本的表达水平,我么你可以利用他们进行下一步的分析
6、差异基因的定义
awk '( $10 <-1 && $12 < 0.05 ){print $0}' gene_exp.diff > up_regulate.txt(上调)
awk '( $10 > 1 && $12 < 0.05 ){print $0}' gene_exp.diff > down_regulate.txt(下调)