RNA-seq流程分析比较之半小时得到差异基因

之前已经写过长篇RNA-seq完整分析过程,本文将介绍采用不同比对软件进行比对，然后采用相同的定量软件进行定量；相同的比对软件，然后采用不同的定量软件进行定量，最后采用目前公认比较靠谱的差异分析R包DESeq2进行差异分析

1、数据下载

• 创建一个包含下载链接的download.txt文件，内容如下

## osdrm2-1
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503154/SRR3503154.sra

## osdrm2-2
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503152/SRR3503152.sra

## DJ-1
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503151/SRR3503151.sra

## DJ-2
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR350/SRR3503149/SRR3503149.sra

- 小贴士 ：ftp链接格式如下
- ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR+SRR号数字的前三位/SRR号/SRR号.sra

2、将sra文件转变为fq文件格式(什么是fq文件格式，先百度，有空再写)

ls *sra |while read id; do nohup fastq-dump --gzip --split-3 $id ; done

3、质量检测

ls *.fastq.gz |xargs nohup fastqc -t 4 -o ./QC &
multiqc ./*zip

4、数据过滤 + 再次质量检测

trim_galore --paired --trim1 -o ./ /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503151_1.fastq.gz /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503151_2.fastq.gz

trim_galore --paired --trim1 -o ./ /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503149_1.fastq.gz /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503149_2.fastq.gz

trim_galore --paired --trim1 -o ./ /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503154_1.fastq.gz /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503154_2.fastq.gz

trim_galore --paired --trim1 -o ./ /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503152_1.fastq.gz /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/row_data_fq/SRR3503152_2.fastq.gz

ls *.fastq.gz |xargs nohup fastqc -t 8 -o ./QC &

- 小贴士：一般公司会给我们一个clean data 和 一个row data，现在测序公司竞争这么激烈，一般是直接可以拿clean data跳过 1-4 步进行下一步分析的。

5、对基因组文件建立索引

bowtie2-bulid all.fa rice  ## bowtie2建立索引

hisat2-build --ss chrX.ss --exon chrX.exon genome/all.fa chrX_tran # hisat2建立索引

subread-buildindex -o rice rice.fa  subread建立索引

6、将我们的测序数据比对到物种参考基因组上

• subjunc

subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice \
-r /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503151_1_val_1.fq.gz \
-R /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503151_2_val_2.fq.gz \
-o DJ-1.bam

subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice \
-r /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_1_val_1.fq.gz \
-R /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_2_val_2.fq.gz \
-o DJ-2.bam

subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice \
-r /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503154_1_val_1.fq.gz \
-R /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503154_2_val_2.fq.gz \
-o drm2-1.bam

subjunc -T 5 -i /public/home/qliu/data/index/rice_subread_index/rice \
-r /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503152_1_val_1.fq.gz \
-R /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503152_2_val_2.fq.gz \
-o drm2-2.bam

• tophat2

tophat2  -p 8 -o DJ-1 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 \
/public/home/qliu/RNA_practice/rice_bowtie2_index/rice \
/public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503151_1_val_1.fq.gz \
/public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_2_val_2.fq.gz

tophat2  -p 8 -o DJ-2 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 \
/public/home/qliu/RNA_practice/rice_bowtie2_index/rice \
/public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_1_val_1.fq.gz \
/public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_2_val_2.fq.gz

tophat2  -p 8 -o drm2-1 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 \
/public/home/qliu/RNA_practice/rice_bowtie2_index/rice \
/public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503154_1_val_1.fq.gz \ /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503154_2_val_2.fq.gz 

tophat2  -p 8 -o drm2-2 -G /public/home/qliu/RNA_practice/rice_gtf/all.gff3 \
/public/home/qliu/RNA_practice/rice_bowtie2_index/rice \ /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503152_1_val_1.fq.gz \
/public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503152_2_val_2.fq.gz

• hisat2

hisat2 -p 4 -x /public/home/qliu/data/index/rice_hisat2_index/rice \
-1 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503151_1_val_1.fq.gz \
-2 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_2_val_2.fq.gz \
-S DJ-1.sam 2>DJ-1_align.log

hisat2 -p 4 -x /public/home/qliu/data/index/rice_hisat2_index/rice \
-1 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_1_val_1.fq.gz \ 
-2 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503149_2_val_2.fq.gz \
-S DJ-2.sam 2>DJ-2_align.log

hisat2 -p 4 -x /public/home/qliu/data/index/rice_hisat2_index/rice \
-1 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503154_1_val_1.fq.gz \ 
-2 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503154_2_val_2.fq.gz 
 \
-S drm2-1.sam 2>drm2-1_align.log

hisat2 -p 4 -x /public/home/qliu/data/index/rice_hisat2_index/rice \
-1 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503152_1_val_1.fq.gz  \ 
-2 /public/home/qliu/data/B512_data/RNA_data/2018_7_23_drm2/trim_data/SRR3503152_2_val_2.fq.gz 
 \
-S drm2-2.sam 2>drm2-2_align.log

samtools sort -n -@ 8 -o HA-1.Nsort.bam HA-1.sam

samtools sort -n -@ 8 -o HA-2.Nsort.bam HA-2.sam

samtools sort -n -@ 5 -o JA-1.Nsort.bam JA-1.sam

samtools sort -n -@ 5 -o JA-2.Nsort.bam JA-2.sam

7、reads的定量分析

- 7-1 对subjunc进行比对的结果采用当前比较热门的定量软件进行分析

• subjunc + featureCounts

featureCounts -p -T 6 -a /public/home/qliu/RNA_practice/rice_gtf/rice.gtf -o DJ_drm2_fCount.out \
/public/home/qliu/drm2-RNA_seq/subjunc_align/DJ-1.bam \
/public/home/qliu/drm2-RNA_seq/subjunc_align/DJ-2.bam \
/public/home/qliu/drm2-RNA_seq/subjunc_align/drm2-1.bam \
/public/home/qliu/drm2-RNA_seq/subjunc_align/drm2-2.bam
## 得到reads的矩阵
awk -F '\t' '{print $1,$7,$8,$9,$10}' OFS='\t' Cr_DJ-osdrm2_fCount.out > WT_osdrm2_matrix.out
第六列为基因的exon长度，可用来计算TPM值用

• subjunc + htseq-count

htseq-count -f bam --stranded=no --type=exon --idattr=gene_id /public/home/qliu/drm2-RNA_seq/subjunc_align/DJ-1.bam /public/home/qliu/drm2-RNA_seq/rice.gtf 1>DJ-1.Counts 2>DJ-1.HTseq.log

htseq-count -f bam --stranded=no --type=exon --idattr=gene_id /public/home/qliu/drm2-RNA_seq/subjunc_align/DJ-2.bam /public/home/qliu/drm2-RNA_seq/rice.gtf 1>DJ-2.Counts 2>DJ-2.HTseq.log

htseq-count -f bam --stranded=no --type=exon --idattr=gene_id /public/home/qliu/drm2-RNA_seq/subjunc_align/drm2-1.bam /public/home/qliu/drm2-RNA_seq/rice.gtf 1>drm2-1.Counts 2>drm2-1.HTseq.log

htseq-count -f bam --stranded=no --type=exon --idattr=gene_id /public/home/qliu/drm2-RNA_seq/subjunc_align/drm2-2.bam /public/home/qliu/drm2-RNA_seq/rice.gtf 1>drm2-2.Counts 2>drm2-2.HTseq.log

## 合并各个Counts文件
paste DJ-1.Counts DJ-2.Counts drm2-1.Counts drm2-2.Counts | awk '{printf $1"\t";for(i=2;i<=NF;i=i+2)printf $i"\t";print $i}' > WT_osdrm2_matrix.out

- 7-2 对不同比对软件比对得到的bam文件采用featureCounts进行定量分析

• tophat2
• hisat2
• subjunc

命令和上面类似，只是将bam文件名改一下而已

8、进行差异分析

• DESeq2

setwd("E://B512/Project/data/drm2_706/my_ananlysis/") ## 设置当前运行路径
countdata <- read.table("WT_osdrm2_matrix.out", header=TRUE, row.names=1) 
colnames(countdata) <- c("DJ-1","DJ-2","drm2-1","drm2-2")
condition <- factor(c("DJ","DJ","drm2","drm2"))
coldata <- data.frame(row.names=colnames(countdata), condition)
dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)
dds <- DESeq(dds)
res <- results(dds)
table(res$padj<0.01)
res_padj <- res[order(res$padj), ]
write.csv(res, file="diffexpr_padj_results.csv",row.names = T)

subset(res,padj < 0.01) -> diff
subset(diff,log2FoldChange < -1) -> down
subset(diff,log2FoldChange > 1) -> up

as.data.frame(down) -> down_gene
as.data.frame(up) -> up_gene
write.csv(up_gene, file="up_gene_LogFC1.csv",row.names = T)
write.csv(down_gene, file="down_gene_LogFC1.csv",row.names = T)

9、结果比较

• 从下图可以看出featureCounts和htseq-count定量软件之间的差别几乎可以忽略，但是前者定量时间短，且不需要合并得到的Counts文件

  
  image.png

• 从下图可以看出这三个比对软件之间的差异几乎很少，但是速度subjunc >>> hisat2 >>..>>> tophat2。

  
  

总结，现在RNA-seq分析工具已经很齐全了，通过今天这几款流程的比较，得出半小时解决RNA-seq从比对到拿到差异基因的最快方法 subjunc + featureCounts + DESeq2。当然细节问题方面，可能没注意。有不得理的地方还望指出。谢谢！