RNAseq分析流程

服务器准备

本实验室有一台服务器,但是由于年代久远(十年前的一个师姐搭建的服务器),导致服务器总是报错,甚至连网络都连不上。因此放弃使用自己实验室的服务器。于是我上淘宝上租了个服务器,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

1683701075275.png

该文件也是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分析的,为了自己的实验进展,我不得不从一个生信小白一步一步的开始学习怎么进行分析。其中一步步的艰辛困苦只有自己知道,有时候抓耳挠腮的不知道怎么解决报错的时候,请教别人,别人一眼就看出来是哪里出现了错误……总之对于我这个小白来说整个过程非常的艰难。好在坚持就是胜利,那些不能击败我的,将使我更强大。在此我把整个流程写下来,方便自己记忆,也给别人提供一个生信小白进行转录组分析的流程

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,185评论 6 503
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,652评论 3 393
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 163,524评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,339评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,387评论 6 391
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,287评论 1 301
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,130评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,985评论 0 275
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,420评论 1 313
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,617评论 3 334
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,779评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,477评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,088评论 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,716评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,857评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,876评论 2 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,700评论 2 354

推荐阅读更多精彩内容