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}.countstep.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.txtPfam注释
使用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.changeKEGG注释
下载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()结语
此次流程梳理没有将结果展示出来，主要是方便快速阅读整理。该流程还有很多地方可以优化，之后会单独列文章进行详细介绍。