2.软件以及流程以及代码: 本实验使用常规的Tophat2 比对
(1)样本制备、建库:总RNA ---> polyA富集mRNA ---> 打断 ---> 随机引物反转录成cDNA ---> 末端修复、加A、加接头
RNA提取:方法-使用trizol-based方法或者试剂盒方法;
富集mRNA,除去rRNA等RNA:依据为-在测mRNA过程中,首先要去除rRNA。以人为例,在抽提的总RNA中,95%的RNA是rRNA,2%的RNA是mRNA,剩下的则是lncRNA、microRNA、siRNA等。rRNA整个人类当中是非常保守的,在各个组织器官中也是非常稳定的,因此这些测序结果对我们的研究是没有用处的。mRNA则是RNA中比较重要的部分。
具体:以人的为例:总RNA跟带有Poly(T)探针的磁珠结合-洗脱结合的mRNA-用Mg离子溶液打碎mRNA-随机引物反转第一条链的cDNA,之后再合成第二条cDNA,获得双链cDNA-对双链cDNA末端修复,加A加接头-片段选择,PCR扩增、纯化(如果样本中存在污染物,则需要结合试剂盒进一步纯化)。
###一般在会测序前对总RNA进行一次质检,根据电泳质检结果中的18S和28S(rRNA)两个峰的高度以及峰的尖度来判断RNA的质量,峰越高越尖(RIN > 8.0)表示RNA的完整度越好。当然,浓度以及A260/A280比值也是需要的。
1、真核生物种常规的去除rRNA的方法是通过oligo(dT)富集带有polyA尾的mRNA来实现的,2、不含有polyA尾的转录本序列以及存在部分降解的总RNA样本,所以这种方法针对福尔马林(Formalin-Fixed)样本和FFPE(Paraffin-Embedded)石蜡包埋样本是不适用的,否则对获得样本中最全面的转录本信息会产生显著影响,一般采用需结合RiboZero、RiboMinus等是结合来开展去除。针对FFPE样本还有结合双链特异性核酸酶构建文库来降低后续测序数据中的rRNA序列比例的。
建库:去除rRNA之后获得的mRNA进行构建文库,先对mRNA打断再进行反转录的文库构建方法。之后反转的cDNA再末端修复到平末端,加上ployA和接头。
###当然,里面涉及的蛋白相关知识,例如蛋白变性失活、键的破坏、复性和盐析等涉及的蛋白四级结构、各种化学键和相互作用的内容。
(2)测序:SE测序-单端测序;PE测序-双端测序;一般现在用的是PE测序。
##里面有测序的方法、测序原理、不同品牌的区别
.sra-> .fastq 代码:fastq-dump --gzip --split-3 –O ../fastq/ -A ../xx.sra
(3)质量控制:当然用IGV也可以做一些质控。
Fastqc 自己命名.fastq.gz -o保存文件夹/
(4)数据与处理:
比对质量过滤、修剪:trimmomatic PE输入文件.1.fastq.gz 输入文件.2.fastq.gz paired1.fq.gz unpaired.1.fq.gz paired2.fq.gz unpaired.2.fq.gz AVGQUAL:20 MINLEN:50(删去质量小于20,且删除读段小于50的片段)
除去那些测不出来或者说未被识别的序列 标记为N,很多的话配对时需要删除这样的序列:prinseq-lite.pl –fastq read1.fastq
–fastq2 read2.fastq –ns_max_n 2 –out_god nfiltered –out_bad null –no_qual_header
– log –verbose (删除每条序列上含有2个N的序列)
去接头:trimmomatic PE输入文件.1.fastq.gz 输入文件.2.fastq.gz paired1.fq.gz unpaired.1.fq.gz paired2.fq.gz unpaired.2.fq.gzILLUMINACLIP:TruSeq2 –PE .fa:2:30:10:1:true (删除TruSeq2接头,允许有2个不匹配,回文剪接阈值是30,简单剪接阈值10,回文模式检测到的最低接头长度是1,反向读段被保留-默认它被删除)
重复:prinseq-lite.pl–fastq read1.fastq –fastq2 read2.fastq –derep_min 101
–out_god nfiltered –out_bad
null –no_qual_header – log –verbose
(5)比对or 的de novo:
tophat2 –o 保存文件夹/ --transcriptome-index 转录本索引bt2/ -p 8 基因组索引bt2/ 输入文件
(6)比对结果注释,量化RSeQc
1.比对统计项:bam_stat.py –i 比对文件
2.比对到基因组各个原件上的情况:
Read_distribution.py –r 基因组的gtf对应的bed文件–I 比对后的文件bam
3.转录本覆盖度是不是有偏差
geneBody-coverage.py –r 基因组的gtf对应的bed文件–I 比对后的文件bam
4.测序深度上表达丰度检测
RPKM_saturation.py –r 基因组的gtf对应的bed文件–I 比对后的文件bam
5. 剪接点
Junction_annotation.py–r 基因组的gtf对应的bed文件–I 比对后的文件bam
6.剪接点饱和状态检测
Junction_saturation.py –r 基因组的gtf对应的bed文件–I 比对后的文件bam
量化:
1. 基因水平上-reads落在那些基因上 HTSeq软件
比对结果按照名字排序:samtools sort –n 比对文件 排序后的文件
量化:htseq-count –f bam –stranded=no 排过序的文件bam基因组索引文件gtf > counts.txt
2. 转录水平量化:cufflinks –G 基因组文件gtf -b基因组–u –p 8 比对后为排序的文件bam -o保存文件夹
3. 外显子水平 DEXSeq软件
扁平化(即,先将基因组注释文件扁平化,拉开距离,形成不重叠的外显子区域,进而将比对数据进行比对)
Python2 dexseq_prepare_annotation.py 基因组索引文件gtf 扁平化文件.gtf
量化 python2 dexseq_count.py –p yes–s no –r name扁平化的文件.gtf 按照名字排序后的sam文件 输出文件.txt
(7)组装:
Cufflinks -P 8 –O保存文件夹/ 比对后的文件bam/
De novo:
Trinity.pl –seqType fq –JM 10G –left 1.fq–-right 2.fq –CPU 4
(8)差异表达分析:R语言作图
参照:http://www.360doc.com/content/18/0309/18/33459258_735717104.shtml http://www.360doc.com/content/17/0911/22/46164085_686360709.shtml http://www.biotrainee.com/thread-1084-1-1.html
setwd("c:/Users/du/Desktop/R/trails/")
##建立数据框 转录组数据比对量化结果
condition1
condition2
rep1
rep2
raw_count
Now if you have some more new interesting idea you will jump it > for example you can usejoin 1.txt 2.txt >1_2.txt Then you will get a mature txt and can be used in frame.OK let try it latter.
raw_count_filter
ENSEMBL
row.names(raw_count_filter)
raw_count_filter
##设计矩阵
group
rownames(group)
names(group)
#标准化矩阵
dds
dds
summary(results(dds))
res
ressult
write.table(ressult,file="gene1_RNAseq2.xls",sep = "\t",quote=F,row.names = T,col.names = NA)
###提取差异分析结果
x
x
attach(x)
x1=1 | log2FoldChange <=-1))
write.table(x1,"differ_out_gene.xls",sep = "\t",quote=F,row.names = T,col.names = NA)
##GO KEGG GSEA analysis ***(clusterProfiles package 可以下载现成的包,再导入。) 还需要DOSE DO.db包
##小鼠的如下: 注释数据在R包中有
library(org.Mm.eg.db)
keytypes(org.Mm.eg.db)
library(clusterProfiler)
library(DOSE)
#GO analysis
#gene_id 转换
gene
transsid
#看单个基因的话
single_Ap1m2
##**正式GO分析
ego
###ego1
summary(ego)
head(summary(ego))
#gorich 后做几个图:
##网络图
emapplot(ego)
#goplot(ego)
goplot(ego)
#Bar plot
barplot(ego, showCategory=20)
# 气泡图
dotplot(ego,font.size=5,showCategory=30)
#Gene-Concept Network
cnetplot(ego)
#UpSet Plot
upsetplot(ego)
#Heatmap-like functional classification
heatplot(ego)
###GSEA分析
#获取按照log2FC大小来排序的基因列表
genelist
names(genelist)
genelist
# GSEA分析(具体参数参考:https://mp.weixin.qq.com/s/p-n5jq5Rx2TqDBStS2nzoQ)
gsemf
head(gsemf)
gseaplot(gsemf,geneSetID = "GO:0000977")
###KEGG(pathway)分析
# 转换ID适合KEGG
m=bitr(rownames(x1),fromType = 'ENSEMBL',toType = 'ENTREZID', OrgDb = 'org.Mm.eg.db')
kegg
# KEGG分析,在KEGG官网中,物种都有对应的缩写,小鼠mmu,其他的缩写看官网:http://www.genome.jp/kegg/catalog/org_list.html
kk
head(summary(kk))
# 气泡图
dotplot(kk, font.size=5)
# 将GO/KEGG结果转换成CSV格式输出
write.table(as.data.frame(kk),"KEGG-enrich.xls",row.names =F)
write.csv(as.data.frame(ego),"GO-enrich.xls",row.names =F)
NOTICE:最近在学习和总结以及做一些自己的课题,因此会时常更新一部分,敬请原谅!
我看了一下:HTSeq软件与cufflinks差异区别在哪里
如下:HTSeq产生的是reads匹配到基因外显子上的序列,得到的是数据库中基因名字(需要在R中转换为俗称可认识的名字,之后获得差异表达基因)图一;cufflinks:如果直接生成的是有FPKM值跟HTSeq类似吧图二;而cufflinks先组装再合并再获得差异基因,之后提取差异基因结果就是免去了R语言中很多的操作图三与图四;如下:
看看R中省去了好多步骤:如下:
后续理解上:
事实上:可以用hisat2代替tophat2,毕竟作者也是这么建议的。
步骤分为如下过程,为了好理解:
1. 数据下载:geo数据库的测序数据;UCSC网站基因组数据(chromFa.tar.gz);gencode网站基因组注释文件(gtf);hisat2网站的index文件;RSeQc软件网站作覆盖度的文件(.bed)
2.数据比对:hisat2 只是用了index文件; samtools的格式转换--排序等
3.比对结果质检:RSeQc的质检;
4.read计数,read归类于哪个基因(HTSeq)、转录(cufflinks)、外显子区域(DEXSeq);
5. 差异表达分析:归类后要进行比较了;DESeq2包
6. 富集分析:因为某一个基因不能仅凭借表达量多少就判断多少了,如果是低表达的突然高表达一点呢,高表达但是比正常状态却少了呢,这时候就需要看看富集到一起的结果是咋样的了。Y叔的clusterfiler包。GO KEGG
7. 其他分析:聚类分析图;主成分分析;(R语言中高级技能包括四类:广义线性模型、聚类分析、时间序列、主成分分析);别的嘛,需要什么就怎么操作吧。