前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >每月一生信流程之rnaseqDTU(差异转录本)

每月一生信流程之rnaseqDTU(差异转录本)

作者头像
生信技能树
发布2019-11-09 17:19:45
2.2K0
发布2019-11-09 17:19:45
举报
文章被收录于专栏:生信技能树

上一期我们推荐的是转录组经典表达量矩阵下游分析大全 本期我们聊聊可变剪切,流程里面写的差异转录本,或者差异外显子,都差不多的意思。

全部bioconductor流程链接在;http://www.bioconductor.org/packages/release/BiocViews.html#___GeneExpressionWorkflow 目前的27个流程,已经分门别类的整理好了,我们每个月学一个流程,预计两年就可以成为生物信息学领域的全栈工程师啦!

一图看懂DGE, DTE and DTU

参考文献:F1000Research 2019, 8:213 Last updated: 18 MAR 2019

image-20191106170904529

今天学习rnaseqDTU

我们首先看看转录组领域的基因表达相关流程吧,首先一起学习的是Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification

  • http://www.bioconductor.org/packages/release/workflows/vignettes/rnaseqDTU/inst/doc/rnaseqDTU.html

这个就没有中文版了,实际上是使用salmon软件把我们的RNA-seq测序的fastq数据根据参考转录组进行定量后,走几个R包流程。

在R里面安装这个bioconductor流程

代码语言:javascript
复制
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("rnaseqDTU")

因为原文写的实在是太详细了,我这里就不拷贝粘贴了,大家直接去阅读即可:

  • RNA-seq workflow for differential transcript usage following Salmon quantification
  • http://www.bioconductor.org/packages/release/workflows/vignettes/rnaseqDTU/inst/doc/rnaseqDTU.R

全部目录如下;

  • 3 Methods
  • 3.1 Simulation
  • 3.2 Operation
  • 3.3 Quantification and data import
  • 3.4 DTU testing
  • 4 Workflow
  • 4.1 Salmon quantification
  • 4.2 Importing counts into R/Bioconductor
  • 4.3 Transcript-to-gene mapping
  • 5 Statistical analysis of differential transcript usage
  • 5.1 DRIMSeq
  • 5.2 stageR following DRIMSeq
  • 5.3 Post-hoc filtering on the standard deviation in proportions
  • 5.4 DEXSeq
  • 5.5 stageR following DEXSeq
  • 5.6 Citing methods in published research
  • 6 DTU analysis complements DGE analysis
  • 6.1 DGE analysis with DESeq2
  • 6.2 DGE analysis with edgeR
  • 6.3 End of workflow section

流程代码

首先是salmon的:

代码语言:javascript
复制
salmon quant -i $index -l A -1 $fq1 -2 $fq2  -p 4 -o quants/${sample}_quant

每个样本的fq测序数据都会被

image-20191106090849429

然后把所有的样本的quant.sf文件批量读取到R里面:

代码语言:javascript
复制
rm(list = ls())  
options(stringsAsFactors = F)
library(rnaseqDTU)
library(tximport)
files=list.files('quants/',pattern = 'quant.sf',recursive = T,full.names = T)
txi <- tximport(files, type="salmon", txOut=TRUE,
                countsFromAbundance="scaledTPM")
cts <- txi$counts
cts <- cts[rowSums(cts) > 0,]

library(GenomicFeatures)
# 文件 gencode.v32.annotation.gtf.gz 自己在gencode数据库网页下载即可
gtf <- "../database/gencode/gencode.v32.annotation.gtf.gz"
txdb.filename <- "gencode.v32.annotation.sqlite"
txdb <- makeTxDbFromGFF(gtf)
saveDb(txdb, txdb.filename)

txdb <- loadDb(txdb.filename)
txdf <- select(txdb, keys(txdb, "GENEID"), "TXNAME", "GENEID")
tab <- table(txdf$GENEID)
txdf$ntx <- tab[match(txdf$GENEID, names(tab))]

cts[1:3,1:3]
range(colSums(cts)/1e6)
head(txdf) 
cts=cts[rownames(cts) %in%  txdf$TXNAME,] 
dim(cts)
txdf <- txdf[match(rownames(cts),txdf$TXNAME),]
all(rownames(cts) == txdf$TXNAME)

counts <- data.frame(gene_id=txdf$GENEID,
                     feature_id=txdf$TXNAME,
                     cts)
save(counts,files,file = 'salmon-out.Rdata')

上面整理salmon结果的代码,看起来很复杂,其实修改的地方不多,值得注意的是:

image-20191106091705316

表达矩阵的第一列是基因的ID,第二列是转录本的ID,后面才是表达量哦。

有一个分组信息,我这里并没有给出我的代码,因为每个人的项目不一样,你需要自己制作,但凡有点R语言基础的,都是没有问题啦。就是 samps 那个变量的内容,有了它,下面的DRIMSeq流程分析差异转录本表达量才有意义。

接着运行DRIMSeq流程即可:

代码语言:javascript
复制
library(DRIMSeq)
d <- dmDSdata(counts=counts, samples=samps)
d
counts(d[1,])[,1:4]
n <- 12
n.small <- 6
d <- dmFilter(d,
              min_samps_feature_expr=n.small, min_feature_expr=10,
              min_samps_feature_prop=n.small, min_feature_prop=0.1,
              min_samps_gene_expr=n, min_gene_expr=10)
d
table(table(counts(d)$gene_id))
design_full <- model.matrix(~condition, data=DRIMSeq::samples(d))
colnames(design_full)
table(samps$condition)
set.seed(1)
system.time({
  d <- dmPrecision(d, design=design_full)
  d <- dmFit(d, design=design_full)
  d <- dmTest(d, coef="conditionControl")
})
res <- DRIMSeq::results(d)
head(res)
res.txp <- DRIMSeq::results(d, level="feature")
head(res.txp)
no.na <- function(x) ifelse(is.na(x), 1, x)
res$pvalue <- no.na(res$pvalue)
res.txp$pvalue <- no.na(res.txp$pvalue)
save(d,res,res.txp,file = 'DRIMSeq-out.Rdata')

是不是非常简单,就拿到了全部的转录本水平的差异表达 呢,还可以可视化如下:

image-20191106091900001

可以看到,我举例的这个项目里面,tumor和control组的样本量是不平衡的,而且基因ID也不容易理解,大家可以自行转换为基因的symbol,这样出图更直观。

学习这样的流程是需要一定背景知识的

首先是LINUX学习

我在《生信分析人员如何系统入门Linux(2019更新版)》把Linux的学习过程分成6个阶段 ,提到过每个阶段都需要至少一天以上的学习:

  • 第1阶段:把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅,主要目的就是去可视化,熟悉黑白命令行界面,可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。
  • 第2阶段:做到文本文件的表格化处理,类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余,查找,切割,替换,合并,补齐,熟练掌握awk,sed,grep这文本处理的三驾马车。
  • 第3阶段:元字符,通配符及shell中的各种扩展,从此linux操作不在神秘!
  • 第4阶段:高级目录管理:软硬链接,绝对路径和相对路径,环境变量
  • 第5阶段:任务提交及批处理,脚本编写解放你的双手
  • 第6阶段:软件安装及conda管理,让linux系统实用性放飞自我
然后是R学习

我在在生信分析人员如何系统入门R(2019更新版) 里面给初学者的知识点路线图如下:

  • 了解常量和变量概念
  • 加减乘除等运算(计算器)
  • 多种数据类型(数值,字符,逻辑,因子)
  • 多种数据结构(向量,矩阵,数组,数据框,列表)
  • 文件读取和写出
  • 简单统计可视化
  • 无限量函数学习
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-11-07,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 一图看懂DGE, DTE and DTU
  • 今天学习rnaseqDTU
  • 在R里面安装这个bioconductor流程
  • 流程代码
  • 学习这样的流程是需要一定背景知识的
    • 首先是LINUX学习
      • 然后是R学习
      领券
      问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档