流程代码

首先是salmon的：

salmon quant -i $index -l A -1 $fq1 -2 $fq2  -p 4 -o quants/${sample}_quant
每个样本的fq测序数据都会被



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

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结果的代码，看起来很复杂，其实修改的地方不多，值得注意的是：





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

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

接着运行DRIMSeq流程即可：

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，这样出图更直观。
※ 修改:·purplesoul 于 Apr  7 15:35:00 2020 修改本文·[FROM: 120.236.174.*]
※ 来源:·水木社区 http://www.newsmth.net·[FROM: 120.236.174.*]

修改:purplesoul FROM 120.236.174.*
FROM 120.236.174.*