如何使用不基于比对的salmon定量流程

测试开头

今天是生信星球陪你的第488天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

豆豆写于19.11.19
电脑没电了，借用花花电脑编辑一下。
今天的内容很简单，就是介绍一下salmon定量软件的上下游使用。

salmon上游分析

它的帮助文档在：https://salmon.readthedocs.io/en/latest/salmon.html

salmon需要对转录本构建索引。因此只能使用参考转录组，而不能使用基因组

构建索引

ref=$wkd/reference/参考转录组.fa# 下载的话可以去GENCODE网站去看看（不过那里只有人和小鼠的）
salmon index -t $ref -i salmon.idx

定量

# 以其中SRR7815790为例
index=$wkd/salmon/salmon.idx
salmon quant -i $index -l A --validateMappings 
        -1 SRR7815790_1.fastq.gz -2 SRR7815790_2.fastq.gz 
        -p 10 -o SRR7815790_quant 1>SRR7815790.salmon.log 2>&1 
# -l 指定文库类型：（A的含义）allow Salmon to automatically infer the library type
# --validateMappings参数的作用：improve both the sensitivity and specificity of mapping and, as a result, can improve quantification accuracy
# -p：线程数

运行速度相比于基于比对的软件（如hisat2）提高很多，例如下面👇

# for hisat2
SRR7815790.hisat.log:Execution time for SRR7815790 hisat and sam2bam : 115.350829 seconds
# for salmon
SRR7815790.salmon.log:Execution time for SRR7815790 salmon : 36.016640 seconds

Salmon的结果

和Hisat的不同，它的每个数据都是独立的文件夹，其中quant.sf 就是每个样本的定量结果

Salmon的下游分析

下面列出了代码，具体的使用重点就是 配置tx2gene 这部分，有了表达矩阵，就能再进行差异分析和后续的操作啦

对于salmon的定量结果，需要用到R里面的 tximport 导入

tximport函数主要需要两个参数：定量文件files和转录本与基因名的对应文件tx2gene

rm(list = ls())
options(stringsAsFactors = F)

####################
# 配置files路径
####################
dir <- file.path(getwd(),'quant/')
dir
files <- list.files(pattern = '*sf',dir,recursive = T)
files <- file.path(dir,files)
all(file.exists(files))

####################
# 配置tx2gene
# 如何构建？以小鼠为例
####################
# https://support.bioconductor.org/p/101156/
# BiocManager::install("EnsDb.Mmusculus.v79")
if(F){
  library(EnsDb.Mmusculus.v79)
  txdf <- transcripts(EnsDb.Mmusculus.v79, return.type="DataFrame")
  mm10_tx2gene <- as.data.frame(txdf[,c("tx_id", "gene_id")])
  head(mm10_tx2gene)
  write.csv(mm10_tx2gene,file = 'mm10_tx2gene.csv')
}
tx2_gene_file <- 'mm10_tx2gene.csv'
tx2gene <- read.csv(tx2_gene_file,row.names = 1)
head(tx2gene)

####################
# 开始整合
####################
library(tximport)
library(readr)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,ignoreTxVersion = T)
names(txi)
head(txi$length)
head(txi$counts)
# 发现目前counts的列名还没有指定

####################
# 添加列名
####################
# 获取导入文件名称的ID，如SRR7815870
library(stringr)
# 先得到SRRxxxx_quant这一部分
n1 <- sapply(strsplit(files,'\/'), function(x)x[12])
# 再得到SRRxxxx这一部分
n2 <- sapply(strsplit(n1,'_'), function(x)x[1])

colnames(txi$counts) <- n2
head(txi$counts)


####################
# 操作新得到的表达矩阵
####################
salmon_expr <- txi$counts
# 表达量取整
salmon_expr <- apply(salmon_expr, 2, as.integer)
head(salmon_expr)
# 添加行名
rownames(salmon_expr) <- rownames(txi$counts)
dim(salmon_expr)

save(salmon_expr,file = 'salmon-aggr.Rdata')