RNA-seq数据分析流程详解 - 卖萌控的博客

如果大家想要了解测序原理，参考文章（测序产生了那么多数据，你知道测序的原理吗？）；如果大家想要了解测序数据格式，参考文章（）

本文介绍RNA-seq的具体分析流程。

1、cutadapt去接头

我们拿到的测序数据一般是带有接头的fastq格式文件，需要用cutadapt把接头去掉。具体代码如下：

#cut NAT sample #-u 20(正值u表示切除R1的前20个碱基) -u -30（负值u表示切除R1的前20个碱基）/ #-U 20(正值U表示切除R2的前20个碱基) -U -30 (负值U表示切除R2的前20个碱基) / #-m 30（丢弃碱基数小于30的reads）-j 6（6个线程，并行加速作用） #-a read1的正向接头序列；-A read2的反向接头序列 nohup cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT / -u 20-u - 30-U 20-U - 30-m 30-j 6/ -o / data2/xxx/ data2/ 2cutadapt/NAT-RNA_combined_clean_R1.fastq / -p / data2/xxx/ data2/ 2cutadapt/NAT-RNA_combined_clean_R2.fastq / / data2/xxx/ data2/rawdata/NAT-RNA_combined_R1.fastq / / data2/xxx/ data2/rawdata/NAT-RNA_combined_R2.fastq >cutadapt2.log 2>& 1&

一般做完去接头会用fastqc查看数据质量是否满足要求。

fastqc结果如下所示：

2、bowtie2建索引

我们需要对下载的参考基因组建立索引，代码如下：

cd /data2/xxx/data2/mm10 nohup bowtie2-build mm10_genome.fa mm10_index >bowtie2-build.log 2>&1 &

mm10的bowtie2结果，第一列为文件大小 3、tophat2比对

把我们的测序序列比对到参考基因组上。代码如下：

#比对NAT sample #-G gtf文件；-o bowtie建立参考基因组索引的公共名（前缀）； nohup tophat2 -p 6-G / data2/xxx/ data2/mm10/gtf/mm10_genes.gtf / -o ./NAT / data2/xxx/ data2/mm10/mm10_index/mm10_genome_bowtie2 / / data2/xxx/ data2/ 2cutadapt/NAT-RNA_combined_clean_R1.fastq / / data2/xxx/ data2/ 2cutadapt/NAT-RNA_combined_clean_R2.fastq >NAT.log 2>& 1&

4、cufflinks定量

对比对后到bam文件定量，即确定每个基因的表达值。代码如下：

#1、cufflinks 转录本组装(NAT sample) #bam文件是tophat2比对输出结果文件 cd / data2/xxx/ data2/ 4cufflinks/NAT nohup cufflinks -o ./ -p 8 ../../ 3tophat/NAT/accepted_hits.bam >NAT.log 2>& 1 #2、cuffmerge 组装后的转录本合并 #assembly_list 是文本文件，里面包含cufflinks的两个gtf结果文件路径 / #/data2/xxx/data2/4cufflinks/CA/trans.gtf / #/data2/xxx/data2/4cufflinks/NAT/trans.gtf cd / data2/xxx/ data2/ 4cufflinks nohup cuffmerge -o ./merge/ assembly_list.txt -p 8> merge.log 2>& 1& #3、cuffquant 转录本定量 #-u 同时对比对上基因组上多个位置的reads进行统计分析 #参考基因组mm10的gtf注释文件；tophat2的比对结果bam文件 cd / data2/xxx/ data2/ 4cufflinks/ 3quantity/NAT nohup cuffquant -p 8-u / data2/xxx/ data2/mm10/gtf/mm10_genes.gtf / / data2/xxx/ data2/ 3tophat/NAT/accepted_hits.bam >NAT.log 2>& 1& #4、cuffdiff 差异分析 #-b 参考基因组mm10的基因组序列 #-u 对多比对位置的read加权然后分配比对到基因组上 #-L 样本名 nohup cuffdiff -b / data2/xxx/ data2/mm10/mm10_genome.fa / / data2/xxx/ data2/ 4cufflinks/ 2merge/merged.gtf / -u -L CA,NAT / data2/xxx/ data2/ 3tophat/CA/accepted_hits.bam / data2/xxx/ data2/ 3tophat/NAT/accepted_hits.bam / -p 8-o ./ >cuffdiff.log 2>& 1&

1、cufflinks组装结果

2、cuffmerge合并结果

3、cuffquant定量结果

4、cuffdiff差异分析结果

5.接下来用r做差异表达基因

##Tools preparation if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2") library("DESeq2") ##The count data preparation mode:character raw_data<-read.table(".../fpkm.txt", header = TRUE) ##convert data.frame to matrix raw_data_mat<-as.matrix(raw_data) ##convert "character" matrix to "numeric" matrix raw_data_mat2<-apply(raw_data_mat, 2, as.numeric) ##delete the first column which is extra raw_data_mat2<-raw_data_mat2[,-1] ##add rownames in the gene name of "raw_data_mat" rownames(raw_data_mat2)<-raw_data_mat[,1] ##check matrix mode is.numeric(raw_data_mat2) ##save "raw_data_mat2" as "count_data" count_data<-raw_data_mat2 ##the number should be integer which the next DESeq2 requires count_data_round<-round(count_data) ##check whether exist duplicated row or it will go wrong anyDuplicated(row.names(count_data) ) ##delete duplicated row in count_data count_data_dedup<-count_data[!duplicated(row.names(count_data)), ] ##check whether exist duplicated row in count_data_dedup anyDuplicated(row.names(count_data_dedup) ) ##make a dds(three replicate) condition<- factor(c("ctrl","ctrl","ctrl", "test","test","test") ) coldata<- data.frame(row.names = colnames(count_data[,1:6]), condition) dds <-DESeqDataSetFromMatrix(round(count_data_dedup), coldata, design= ~ condition) ##DESeq dds2<-DESeq(dds) #resultsNames(dds2) res<-results(dds2) ##pick differential genes or use " diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1) " diff_gene_deseq2<-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) ##merge differential gene data and normalized raw_count data result_data<-merge(as.data.frame(diff_gene_deseq2), as.data.frame(counts(dds2, normalized=TRUE)), by="row.names", sort=TRUE) ##save result file write.csv(result_data,file = "..../ctrl_test")

如果大家不仅想要差异表达基因，还想知道样本的SNP和INDEL,可参考这篇文章（）

拿到SNP和INDEL后如果做注释，可参考这篇文章