DESeq做转录组分析的代码 转录组分析流程

RNA-Seq生信分析全流程

• 摘要
• 第一部分

• step.1 下载数据
• step.2 数据质控

• 第二部分

• step.3序列比对
• step.4 计算基因表达量
• step.5 插入片段长度检验
• step.6 基因表达量从count值转换为FPKM值

• 使用基因组注释，通过R工具包GenomicFeatures获得exon length
• 求reads 总数

• 第三部分

• step.7 进行各样品分析

• 样品间相关性分析
• 各样品FPKM箱线图
• 各样品FPKM密度分布对比图

• step.8 差异表达分析

• 有生物学重复*
• 无生物学重复
• 绘制火山图（ggplot2）

• step.9 差异基因功能注释

• 获取差异基因注释信息

• 比对基因组获得基因名称
• 根据位置信息获取基因序列
• 使用bedtools getfasta生成差异基因序列diff.fa
• 对各个数据库进行功能注释

• swiss_port注释

• Pfam注释

• KEGG注释

• GO注释
• eggNOG注释

• COG注释

• 结语

摘要

上周崴脚，一周在家办公，效率较低，没有输出文章，今天来个大的。本篇文章为入职后第一次完成RNA-Seq全流程项目分析后进行的内容梳理，有一些地方依靠师姐（公司导师）已经完成的perl脚本，个人并不是很能看懂，部分细节还没有完全打通。之后可能会自己用python（也可能是shell）重写之后再慢慢补充。

第一部分

step.1 下载数据

百度网盘下载数据
查看gz压缩格式的内容：zcat 18134R-49-01_S0_L001_R1_001.fastq.gz z | head -n 10

step.2 数据质控

批量去除接口

mkdir trim
for i in *_R1_001.fastq.gz;
do
i=${i%_R1_001.fastq.gz*};
trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o ./trim ${i}_R1_001.fastq.gz ${i}_R2_001.fastq.gz
done

生成质控图片

for i in *_001_val_1.fq.gz;
do
i=${i%_001_val_1.fq.gz*};
fastqc ${i}_001_val_1.fq.gz
unzip ${i}_001_val_1_fastqc.zip
mv ${i}_001_val_1_fastqc/Images/per_base_quality.png ${i}_per_base_quality.png
mv ${i}_001_val_1_fastqc/Images/per_base_sequence_content.png ${i}_per_base_sequence_content.png
rm -rf ${i}_001_val_1_fastqc
done

第二部分

step.3序列比对

与参考基因组比对后，生成sam文件

for i in *_001_val_1.fq.gz;
do
i=${i%_001_val_1.fq.gz*};
hisat2-2.0.4/hisat2 -p 8 -x /home/ref/gencode/ucscindex/hg38/genome -1 ${i}_1_val_1.fq.gz -2 ${i}_2_val_2.fq.gz -S ../../02.align/${i}_align.sam 2> ../../02.align/${i}_align.log
cd ../../02.align

对sam文件进行排序，转换成bam文件并建立索引
索引文件格式为bam.bai

samtools view -bS ${i}_align.sam > ${i}.tmp.bam
rm -rf ${i}_align.sam
samtools sort ${i}.tmp.bam -o ${i}.bam
rm -rf ${i}.tmp.nam
samtools index ${i}.bam
done

比对后进行比对结果需要自己写脚本进行统计

step.4 计算基因表达量

htseq计算表达量
此处只计算单个样品（bam文件）count值，如果要合并count还需要自己写脚本合并

mkdir htseq
for i in `ls *.bam`
do
i=${i%.bam}
htseq-count --stranded=no -f bam ${i}.bam /home/xx/RNA-seq_ref/04_36_homo_2019//ref/gencode/gencode.v29.annotation.gtf > htseq/${i}.count

step.5 插入片段长度检验

使用picard CollectInsertSizeMetrics可以获得插入片段长度的txt文档及PDF图片
运行代码

java -jar picard.jar CollectInsertSizeMetrics \
I=input.bam \
O=insert_size_metrics.txt \
H=insert_size_histogram.pdf \

step.6 基因表达量从count值转换为FPKM值

强调：获得的FPKM文本需要转换换行符，不能使用windows格式下的！！！一定要使用Unix格式，win10笔记本或者notepad在右下角都会显示其文本格式。
获得参数

使用基因组注释，通过R工具包GenomicFeatures获得exon length

运行代码

#/home/songyumo/RNA-seq_ref/04_36_homo_2019//ref/gencode/gencode.v29.annotation.gtf
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v29.annotation.gtf",format="gtf")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
n=t(as.data.frame(exons_gene_lens))
write.

求reads 总数

去除第一行后读取总共行数
根据count与FPKM转换公式编写转换脚本，这里有详细介绍 公式
运算代码（R语言）

#根据参考基因组获得每条reads的exon length /home/ref/gencode/gencode.v29.annotation.gtf
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("/home/ref/gencode/gencode.v29.annotation.gtf",format="gtf")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
n=t(as.data.frame(exons_gene_lens))
#导入所有reads的count值，header = T保证后续计算不会将列名算入而造成错误。如果file和脚本不在同一个路径，记得添加file的绝对路径。 
count <- read.table(file="all_count.txt",header = T) #从第二列开始，将每个样本的count循环转化成FPKM 
i <- 2 
repeat{ 
mapped_reads <- sum(count[1:58721,i])#计算每个样本的mapped reads数,由于count文本最后几行是综述数据，需要舍去，因此选择58721行之前的基因数据。 
FPKM <- count[1:58721,i]/(10^-9*mapped_reads*n[1:58721,1])#计算FPKM值，10^-9是单位转换后的换算，不能放在分母外的原因是mapped_reads*n[1:58721,1]的乘积很大会报错，因此需要先缩小倍数再进行计算 
count = data.frame(count[1:58721,],FPKM)#添加到count表格后面i 
i <- i + 1 
if(i >41){ 
break 
}#样品总共有40个样品，我们需要获取从第2列到第41列count值 
} 
#生成表格列名称 
count_colname <- read.table("all_count.txt",header = F,nrow = 1,as.is=TRUE) 
FPKM_colname <- paste(count_colname[1,],"_FPKM",sep="") 
colname <- c(count_colname,FPKM_colname) 
col_name <- colname[-which(colname=="gene_FPKM")]#删除gene_FPKM 
names(count) <- col_name #生成表格 
write.csv(count,"ALL_FPKM.csv",row.names = FALSE)#row.names是为了去除第一列自动生成的行名，同理col.names=FALSE可以去除第一行自动生成的列名 head(read.csv("ALL_FPKM.csv"))

第三部分

step.7 进行各样品分析

样品间相关性分析

使用R包pheatmap
运行代码

library(pheatmap)
#data=read.csv("generead_counts.csv",header = T,row.names = 1)
data=read.table("all_count.txt",header = T,row.names = 1)
matrix=cor(data)
write.csv(matrix,"gene_coefficient_matrix.csv")
pdf(file="gene_cor.pdf")
pheatmap(matrix,cluster_rows=T,cluster_cols = T,display_numbers = T,fontsize = 5)
dev.off()

各样品FPKM箱线图

运行脚本

library(RSkittleBrewer)
library(genefilter)
library(dplyr)
library(devtools)
pheno_data = read.table("sample.csv",header=T,sep=",",colClasses=c("character","factor"))
fpkm = read.table("FPKM_table.csv",header=T,sep=",")
#fpkm = read.table("genes.fpkm_table",header=T,sep="\t") 
fpkm = fpkm[,c("AbiLN.TFD","AbiLN.PBS")]
#row.names(fpkm)=fpkm[,1]
#fpkm = fpkm[,-1]
fpkm<-fpkm[which(rowSums(fpkm) > 0),]
fpkm = log10(fpkm)
#fpkm = log2(fpkm)
pdf("box1.pdf")
boxplot(fpkm,las=2,ylab='log10(FPKM)', col=rainbow(12),cex.axis=0.5)
#boxplot(fpkm,col=as.numeric(pheno_data$treat1),las=2,ylab='log2(FPKM+0.00001)')
dev.off()

各样品FPKM密度分布对比图

使用R包ggplot2
运行代码

library(ggplot2)
library(plotly)
data<-read.table("FPKM_table.txt",header = T,row.names = 1)
pdf("density.pdf")
#png("density.png")
data<-data[which(rowSums(data) > 0),]
data<-log10(data)
data<-stack(data)
#p1<-ggplot(data,aes(x=data$))
p<-ggplot(data, aes(x = values)) +
scale_colour_hue("Sample")+
theme_bw()+
theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+
#theme(legend.title=element_blank()) +
geom_density(aes(group = ind, color = ind, fill=ind),alpha=0.05)+
guides(fill=FALSE)+
labs(title = "Gene expression density", x="log10(FPKM)", y="density")
#stat_density(aes(group = ind, color = ind),position="identity",geom="line")
#ggplotly(p)
p
dev.off()

step.8 差异表达分析

有生物学重复*

使用R包Deseq2，因为本次项目无生物学重复，这里不详细介绍

无生物学重复

使用R包edgeR进行差异表达分析
运行代码

library("edgeR")
library("limma")
rm(list=ls())
countsall<-read.table("all_count.txt", header=TRUE, row.names = 1, sep="\t",encoding ='UTF-8',check.names = FALSE )
counts <- countsall[,c("AbiLN.TFD","AbiLN.PBS")]
group<-factor(c("AbiLN.TFD","AbiLN.PBS"))
y1 <- DGEList(counts=counts,genes = rownames(counts), group=group)
keep <- rowSums(cpm(y1)>1) >= 2
y1 <- y1[keep,]
y<-calcNormFactors(y1)
count_r <- cpm(y,normalized.lib.sizes=TRUE)
write.csv(count_r,file = "generead_counts.csv")
bcv = 0.1
#bcv = 0.4
#bcv = 0.01
et <- exactTest(y, dispersion=bcv^2)
results = et$table
results$category=counts$cateloge
results$QValue = p.adjust(results$PValue, method = 'fdr')
resall <- merge(as.data.frame(results), as.data.frame(count_r),by="row.names",sort=FALSE)
write.csv(resall,"geneDE_results.edger.csv")
signresdata<-subset(resall,resall$QValue<0.05 & abs(resall$logFC) > 1 )
write.csv(signresdata,file = "geneDEsign_results.csv")

绘制火山图（ggplot2）

运行代码

library(ggplot2)
rm(list=ls())
resdata <- read.csv("geneDE_results.edger.csv",header = T,row.names = 1)
resdata <- na.omit(resdata)
threshold <- as.factor(ifelse(resdata$QValue < 0.05 &abs(resdata$logFC) >= 1 ,ifelse(resdata$logFC >= 1 ,'Up','Down'),'Not'))
pdf(file="geneVolcano.pdf")
ggplot(resdata,aes(x=logFC,y=-log10(QValue),colour=threshold)) +
xlab("log2(Fold Change)")+ylab("-log10(QValue)") +
geom_point(size = 0.5,alpha=1) +
ylim(0,200) + xlim(-8,8) +
scale_color_manual(values=c("green","grey", "red"))
dev.off()

火山图出图效果还不是很好看，需要再继续优化

step.9 差异基因功能注释

获取差异基因注释信息

比对基因组获得基因名称

与人的参考基因组信息（gencode.gtf）进行比对索引,获取注释diff.gtf
运行代码

perl /home/diff2gtf.pl ../geneDEsign_results.csv /home/ref/gencode/chroms/gencode.gtf difftmp.gtf
awk '( ~ /gene/){print}' difftmp.gtf > diff.gtf
rm -rf difftmp.gtf
awk -F "\"" '{print }' diff.gtf > diffid.list

根据位置信息获取基因序列

获取差异基因位置信息
运行代码

perl -e '{my $diff=shift;open IN, "$diff";my %hash;while (){chomp; $hash{$_}=1;}close IN; my $in=shift;open IN,"$in";while (my $a=){chomp$a; my $gene=(split/\t/,$a)[3];if ($hash{$gene}){print "$a\n";}}close IN;}' diffid.list /home/ref/gencode/anno/geneID.bed > diffID.bed

使用bedtools getfasta生成差异基因序列diff.fa

运行代码

bedtools getfasta -name -fi /home/ref/gencode/chroms/chr.fa -bed diffID.bed -fo diff.fa

数据库功能注释
下载人类全基因组注释文本，根据差异基因ID与文本比对进行注释,这里的perl脚本是将count、FPKM值与注释合并在一起，可以自己写脚本合并数据列。
运行代码

perl /home/anno.pl ../geneDEsign_results.csv /home/All_Database_annotation.txt count.txt genes.fpkm_table anno.DEG.txt

对各个数据库进行功能注释

swiss_port注释

使用blastx和swiss本地工具进行比对，这个还不是太熟悉。
运行代码

blastx -d swissprot -q ../diff.fa -k 1 -e 0.00001 -o swiss_dia_matches.m8
perl /home/filtertab.pl swiss_dia_matches.m8 > filterswiss_dia_matches.m8
sed -i "1iquery acc.ver\tsubject acc.ver\t% identity\talignment length\tmismatches\tgap opens\tq. start\tq. end\ts. start\ts. end\tevalue\tbit score" filterswiss_dia_matches.m8
perl /home/anno/swiss_prot/add_anno.pl /home/database/swissprot/info.txt filterswiss_dia_matches.m8 swiss.m8 > swiss_anno.txt

Pfam注释

使用Genebank工具gmhmme2进行处理
运行代码

perl /home/anno/Pfam/getfa.pl ../diff.fa ./
cd fasta
for i in *.fa
do
i=${i%.fa};
#echo ${i}
gmhmme2 -m /root/Software/genemark_hmm_euk_linux_64/ehmm/mtx/human_49_99.mtx -o gms.out.faa.${i} ${i}.fa
done
perl /home/anno/Pfam/gene2protein.pl fasta/ > gms.out.faa
##!!!!!!! vi s/gene /gene_/g
hmmscan -o out.txt --tblout out.tbl --noali -E 1e-5 /home/database/Pfam/Pfam-A.hmm gms.out.faa
grep ">" gms.out.faa > gms2fa.txt
###s/>//g
perl /home/gms2gene.pl gms2fa.txt out.tbl > out.tbl.change

KEGG注释

下载KO.list,以及.keg格式文件
KEGG Automatic Annotation Server (genome.jp)

GO注释

绘制GO分析图
WEGO在线分析工具：GO (genomics.org.cn) 这个在线工具生成的图片只有两种颜色，BP,MF,CC之间没有颜色区分，看看能不能再优化一下

eggNOG注释

进行eggNOG注释
在线分析工具：eggNOG-mapper (embl.de)

COG注释

使用ggplot2绘制COG图像
运行代码

args<-commandArgs(T)
if(length(args)<2){
cat("[usage:] ")        quit("no")}
cazyclass<-args[1]
outfile<-args[2]
library(ggplot2)
dat<-read.table(cazyclass,header=T,sep='\t')
pdf(file=paste(outfile,".pdf",sep=""),width=6*1.5,height=6)
colnames(dat)<-c("Motif","Description","Number")
ggplot(dat,aes(x=dat$Motif,y=dat$Number,fill=paste(dat$Motif,':',dat$Description)))+geom_bar(stat="identity")+theme(legend.title=element_blank(),legend.text=element_text(size=6),legend.key=element_rect(size=6))+labs(x='Function class',y='Number of matched genes')
dev.off()

结语

此次流程梳理没有将结果展示出来，主要是方便快速阅读整理。该流程还有很多地方可以优化，之后会单独列文章进行详细介绍。