c语言处理rna序列,RNAseq 完整操作流程以及后续例子操作

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语言中高级技能包括四类：广义线性模型、聚类分析、时间序列、主成分分析)；别的嘛，需要什么就怎么操作吧。