服务器准备
本实验室有一台服务器,但是由于年代久远(十年前的一个师姐搭建的服务器),导致服务器总是报错,甚至连网络都连不上。因此放弃使用自己实验室的服务器。于是我上淘宝上租了个服务器,139元大宽带美国服务器,八核处理器处理数据也很快。
分析用软件的安装
服务器中安装的软件:
质量控制软件 : fastp
比对软件 : hisat2
Bam 文件处理软件 : samtools
转录组拼接软件: stringtie
以上这些软件都可以用conda安装,conda的安装教程(https://blog.csdn.net/weixin_64316191/article/details/127435729)
另外安装stringtie软件的时候总是报错,网上搜索原因是因为需要虚拟环境(我也不知道是什么意思),根据教程把虚拟环境搭建好了(https://blog.csdn.net/qq_43391414/article/details/124398631)以后每次分析数据我都会先切换到这个虚拟环境再进行数据分析。
自己电脑上需要安装的软件:
电脑服务器文件传送软件 WinSCP
服务器连接软件 PuTTY(这里面的复制粘贴分别是鼠标左键和右键)
差异分析软件 R和Rstudio
文件准备
全基因组序列文件
全基因组注释文件
这两个文件可以从Phytozome网站上下载,分别以 fa 和 gff3 后缀为结尾
公司测序返回的自己的实验数据结果,文件很大,而服务器的存储很小,只有100G,因此我选择每个样品分开分析,这样可以分析一个后删掉上一个没用的文件,不会出现因为硬盘不够而报错的问题。
分析流程
第一步: 数据质量控制
公司返回的数据一般是pair reads 一个样本里有两个文件。
fastp -i HDZ_MOCK_1_1.fq.gz -o HDZ_MOCK_1_1_clean.fq.gz -I HDZ_MOCK_1_2.fq.gz -O HDZ_MOCK_1_2_clean.fq.gz -q 3 -u 50 --length_required 100 -j MOCK_1_1.json -h MOCK_1_1.html
# sameple_name换成自己样品的名字,一定与文件名字一样,名称中间的分隔最好用下划线_
# 对于生成的html 文件,建议查阅里面的内容,从数据被过滤掉的比例,和插入片段的长度分布两个角度对文库的质量有个初步评估
第二步:参考基因组索引(index)的构建
hisat2-build Sitalica_312_v2.fa Sitalica
#hisat2-build 全基因组序列文件 序列文件别名
#这个序列文件需要跟自己的数据在同一个文件夹内
第三步:样品reads 比对到基因组
hisat2 -p 20 -x Sitalica -1 HDZ_MOCK_1_1_clean.fq.gz -2 HDZ_MOCK_1_2_clean.fq.gz -S HDZ_MOCK_1_rep1.sam
# Sitalica是我序列文件的别名, HDZ_MOCK_1_1_clean.fq.gz和HDZ_MOCK_1_2_clean分别是fastp分析后的文件
第四步:映射到基因组中的reads的排序及质量的过滤
生成的sam文件非常大,大概30多GB
samtools sort -@ 8 -o HDZ_MOCK_1_rep1.sorted.bam HDZ_MOCK_1_rep1.sam
# @是线程数,一般几核处理器就用几个线程数
samtools view -bS -h -q 20 -@ 8 HDZ_MOCK_1_rep1.sorted.bam > HDZ_MOCK_1_rep1.q20.sorted.bam
samtools index HDZ_MOCK_1_rep1.q20.sorted.bam
第五步:mapping reads的统计分析
samtools flagstat -@ 5 HDZ_MOCK_1_rep1.sorted.bam > HDZ_MOCK_1_rep1.sorted.bam.txt
samtools flagstat -@ 5 HDZ_MOCK_1_rep1.q20.sorted.bam > HDZ_MOCK_1_rep1.q20.sorted.bam.txt
第六步:基因最长转录本的FPKM计算及提取
stringtie -e -B -p 20 -G Sitalica_312_v2.2.gene.gff3 -o HDZ_MOCK_1_rep1.stringtie.gtf -A HDZ_MOCK_1_rep1.tab HDZ_MOCK_1_rep1.q20.sorted.bam
每一个数据都重复一遍上面的流程,直到所有的文件都分析完成。
第七步:使用DESeq2进行差异表达计算
需要准备两个文件
准备文件1:WT_MOCK_VS_HDZ_MOCK_sample_lst.txt
WT_MOCK_1_rep1 WT_MOCK_1_rep1.stringtie.gtf
WT_MOCK_2_rep1 WT_MOCK_2_rep1.stringtie.gtf
HDZ_MOCK_1_rep1 HDZ_MOCK_1_rep1.stringtie.gtf
HDZ_MOCK_2_rep1 HDZ_MOCK_2_rep1.stringtie.gtf
Note:该文件需要在csv文件中填写每一个名称,再新建txt文件夹,将填好的名称复制到txt中。否则接下来程序无法识别这个文件
准备文件2:prepDE.py
该文件是DESeq官网上下载的文件
官网网址(http://ccb.jhu.edu/software/stringtie/)
下载连接(http://ccb.jhu.edu/software/stringtie/dl/)找到prepDE.py以后
右键另存为就下载好了
第八步:分析差异基因(此后在自己电脑上用Rstudio分析)
需要准备文件:HDZ_PEG_VS_HDZ_MOCK_samples.txt
该文件也是csv文件写好后再复制粘贴到txt文件中,其中condition在第二列的最上面。
#需要在Rstudio上安装DESeq2包
library('DESeq2')
setwd("G:/R_analysis") #这是文件读取的目录
countData <- as.matrix(read.csv("HDZ_PEG_VS_HDZ_MOCK_gene_count_matrix.csv", row.names="gene_id"))
colData <- read.csv("HDZ_PEG_VS_HDZ_MOCK_samples.txt", sep="\t", row.names=1)
dds <- DESeqDataSetFromMatrix(countData, colData, design = ~condition)
dds<-DESeq(dds)
dds
resultsNames(dds)
#使用deseq2包中的result()函数提取deseq2差异分析结果
res=results(dds,contrast = c("condition","HDZ_PEG","HDZ_MOCK"))
head(res)
summary(res)
write.csv(res,file = "HDZ_PEG_VS_HDZ_MOCK_results.csv")
#分析结果可视化与导出
table(res$padj<0.05)
plotMA(res,ylim=c(-2,2))
mcols(res,use.names = TRUE)
#提取差异基因
res <- res[order(res$padj),]
resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
deseq_res<-data.frame(resdata)
up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上调差异表达基因
down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下调差异表达基因
write.csv(up_diff_result,"up_HDZ_PEG_VS_HDZ_MOCK_diff_results.csv") #输出上调基因
write.csv(down_diff_result,"down_HDZ_PEG_VS_HDZ_MOCK_results.csv") #输出下调基因
plot(res$log2FoldChange,res$padj)#绘制火山图
library("pheatmap")
select<-order(rowMeans(counts(dds,normalized=TRUE)),decreasing = TRUE)[1:100]
colData(dds)
df<-as.data.frame(colData(dds)[,c("condition","sizeFactor")])
mycolors=colorRampPalette(c("orange","white","purple"))(100)
pheatmap(assay(vsd)[select,],cluster_rows=FALSE,show_rownames=FALSE,border_color=NA,color=mycolors,cluster_cols=FALSE,annotation_col=df)
library(ggplot2)
volcano<-subset(dealdataNA,select = c(Pvalue,Foldchange))
threshold<-as.factor((log2(volcano$Foldchange)>1.5|log2(volcano$Foldchange)<(-1.5))&volcano$Pvalue<0.05)
r03=ggplot(volcano,aes(log2(Foldchange),-log2(Pvalue),colour=threshold))+geom_point()+xlim(-10,10)
r04=r03+geom_vline(xintercept=c(-1.5,1.5),linetype="dotted",size=1)+geom_hline(yintercept=-log2(0.05),col="blue")
r05=r04+labs(title="Volcanoplot")+theme(plot.title = element_text(hjust = 0.5))
第九步:进行差异基因的GO分析
setwd("G:/GO_analysis")
library(GO.db)
godb <- select(GO.db, keys(GO.db), columns(GO.db))
write.csv(godb, "go2name.csv",row.names=F)
library("clusterProfiler")
# 导入基因列表
gene <- read.csv("down_WT_PEG_VS_HDZ_PEG_diff_GO.csv",header=T)
head(gene)
gene <- as.factor(gene$GeneID)
# 导入注释文件
term2gene <- read.csv("Setaria italica基因组GO注释文件.csv",header=F,sep=",")
term2name <- read.csv("go2name.csv",header = F,sep=",")
# 富集分析
x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.05)
# 输出结果
write.csv(x,"down_WT_PEG_VS_HDZ_PEG_diff_GO.csv",row.names=F)
# 绘制条形图
barplot(x)
# 绘制气泡图
dotplot(x)
该方法参考(https://www.jianshu.com/p/4759cd1c751a)
其中准备的三个文件中第一个待富集文件的第一列一定是GeneID。否则会报错
写在最后
我们实验室导师要求我们自己分析转录组的数据,因为这样的话每个样品能够节省一百块钱。但是我们实验室没有一个师兄师姐会RNAseq分析的,为了自己的实验进展,我不得不从一个生信小白一步一步的开始学习怎么进行分析。其中一步步的艰辛困苦只有自己知道,有时候抓耳挠腮的不知道怎么解决报错的时候,请教别人,别人一眼就看出来是哪里出现了错误……总之对于我这个小白来说整个过程非常的艰难。好在坚持就是胜利,那些不能击败我的,将使我更强大。在此我把整个流程写下来,方便自己记忆,也给别人提供一个生信小白进行转录组分析的流程