首先要知道alevin屬于salmon的一個(gè)工具饲宛，用于定量單細(xì)胞轉(zhuǎn)錄本皆愉，使用的比對(duì)方式是類(lèi)似于kallisto的pseudo align，即不將reads比對(duì)到基因組上艇抠，該算法著重于確定一個(gè) read 屬于哪一個(gè)基因幕庐，而不關(guān)心這個(gè) read 在基因上的位置。（Kallisto: 一個(gè)RNA-seq數(shù)據(jù)快速量化軟件 http://blog.sciencenet.cn/blog-656335-984247.html家淤；Kallisto ‘Pseudoalignment’ 原理：https://tinyheero.github.io/2015/09/02/pseudoalignments-kallisto.html）
因此最后沒(méi)法生存單細(xì)胞比對(duì)的bam文件异剥，這個(gè)bam文件還是有很多用處的，比如看一看單細(xì)胞中的突變位點(diǎn)情況絮重。雖然alevin用下來(lái)確實(shí)很快冤寿，但是歹苦，不符合我的分析要求~
alevin的參考文檔：
https://salmon.readthedocs.io/en/latest/alevin.html
使用Salmon對(duì)單細(xì)胞轉(zhuǎn)錄本進(jìn)行定量
下載必要的文件
#### 僅包括蛋白編碼的轉(zhuǎn)錄本
wget -nv ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/gencode.v35.pc_transcripts.fa.gz
#### 包括全部轉(zhuǎn)錄本
wget -nv ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/gencode.v35.transcripts.fa.gz
下載基因組注釋gtf
wget -nv ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/gencode.v35.primary_assembly.annotation.gtf.gz
生成ENST到ENSG的匹配
bioawk -c gff '$feature=="transcript" {print $group}' <(gunzip -c gencode.v35.primary_assembly.annotation.gtf.gz) | awk -F ' ' '{print substr($4,2,length($4)-3) "\t" substr($2,2,length($2)-3)}' - > txp2gene.tsv
生成ENST到Gene Symbol的匹配
bioawk -c gff '$feature=="transcript" {print $group}' <(gunzip -c gencode.v35.primary_assembly.annotation.gtf.gz) | awk -F ' ' '{print substr($4,2,length($4)-3) "\t" substr($8,2,length($8)-3)}' > txp2gene_symbol.tsv
生成salmon的索引序列
./salmon-latest_linux_x86_64/bin/salmon index \
-i transcript_index \
-k 31 \
--gencode \
-p 4 \
-t gencode.v35.transcripts.fa.gz
運(yùn)行單細(xì)胞轉(zhuǎn)錄本定量工具Alevin
./salmon-latest_linux_x86_64/bin/salmon alevin -l ISR \
-1 Sample_L001_R1_001.fastq.gz Sample_L002_R1_001.fastq.gz Sample_L003_R1_001.fastq.gz Sample_L004_R1_001.fastq.gz Sample_L005_R1_001.fastq.gz Sample_L006_R1_001.fastq.gz Sample_L007_R1_001.fastq.gz \
-2 Sample_L001_R2_001.fastq.gz Sample_L002_R2_001.fastq.gz Sample_L003_R2_001.fastq.gz Sample_L004_R2_001.fastq.gz Sample_L005_R2_001.fastq.gz Sample_L006_R2_001.fastq.gz Sample_L007_R2_001.fastq.gz \
--dropseq \
-i transcript_index \
-p 10 \
-o transcripts_output \
--tgMap txp2gene_symbol.tsv \
--expectCells 8000
速度確實(shí)很快，dropseq alignment protocol需要1天才能完成督怜，alevin花了1個(gè)小時(shí)殴瘦。
但是定量結(jié)果的差異還是很大：
dropseq alignment protocol最后比對(duì)出了8000個(gè)左右的細(xì)胞叉钥，
alevin比對(duì)出4300個(gè)細(xì)胞（其中還包括1000個(gè) low confidence）敛惊，這個(gè)原因下文會(huì)初步探索一下。
dropseq alignment protocol利用STAR 2-pass mode比對(duì)率在80%荒澡，
alevin比對(duì)率45%左右姨蟋，這一點(diǎn)在Alevin github issue中有人提問(wèn)過(guò)屉凯，開(kāi)發(fā)者給出的回答是：STAR（或Hisat2）給出的比對(duì)率是全基因組的比對(duì)率，比對(duì)上的reads不一定是轉(zhuǎn)錄本芬探，而alevin的比對(duì)率是reads比對(duì)至基因組轉(zhuǎn)錄本的比率神得。
差異可能來(lái)自于：
- 轉(zhuǎn)錄本的比對(duì)和技術(shù)策略：比對(duì)和“假比對(duì)”的技術(shù)差異
- alevin對(duì)barcode和umi有比較好的校正算法，drop-seq protocol的校正比較簡(jiǎn)單
- 選擇的參考序列和注釋文件的差異：drop-seq protocol使用全基因組序列建立索引偷仿，而alevin僅使用轉(zhuǎn)錄本序列
- 其他原因
后續(xù)的探索
比較一下alevin和drop-seq protocol的比對(duì)結(jié)果
alevin_result <- data.table::fread(input = "featureDump_alevin_result_20200929.txt")
# 讀入drop-seq protocol計(jì)數(shù)的矩陣
p <- data.table::fread("./CH4-LN_hg38_gene_exon_tagged_CODING_UTR_INTRONIC.dge.txt.gz", data.table = F)
# alevin計(jì)數(shù)到的細(xì)胞中有98%是在drop-seq protocol中的
sum(alevin_result$CB %in% colnames(p))/nrow(alevin_result);rm(p)
# 取兩次比對(duì)結(jié)果的交集
p <- FetchData(All, c("nCount_RNA","nFeature_RNA")) %>%
tibble::rownames_to_column("CB") %>%
inner_join(.,alevin_result, by = "CB")
library(ggplot2)
library(grid)
library(gridExtra)
library(ggpubr)
# 比較兩個(gè)數(shù)據(jù)集reads數(shù)的比對(duì)結(jié)果
fig1 <- ggplot(p, aes(x=nCount_RNA, y=DeduplicatedReads)) +
geom_point() +
ggtitle("Dropseq-Tools VS Alevin") +
geom_smooth(method=lm, se=FALSE) +
scale_x_continuous(name = "Dropseq-Tools", limits = c(1E3,3E4)) + #, breaks = seq(5, 15, 2)
scale_y_continuous(name = "Alevin", limits = c(1E3,3E4)) + # , breaks = seq(5, 15, 2)
# annotation_custom(grob1) +
theme(plot.title = element_text(hjust = 0.5)) +
stat_cor(method="pearson") +
theme_light()
fig1
# 比較兩個(gè)數(shù)據(jù)集基因數(shù)的比對(duì)結(jié)果
fig2 <- ggplot(p, aes(x=nFeature_RNA, y=NumGenesExpressed)) +
geom_point() +
ggtitle("Dropseq-Tools VS Alevin") +
geom_smooth(method=lm, se=FALSE) +
scale_x_continuous(name = "Dropseq-Tools", limits = c(300,5E3)) + #, breaks = seq(5, 15, 2)
scale_y_continuous(name = "Alevin", limits = c(300,5E3)) + # , breaks = seq(5, 15, 2)
# annotation_custom(grob1) +
theme(plot.title = element_text(hjust = 0.5)) +
stat_cor(method="pearson") +
theme_light()
fig2
# 看在Alevin中比對(duì)得到的細(xì)胞在Dropseq-Tools數(shù)據(jù)中的分布情況
All$inAlevin <- (Cells(All) %in% alevin_result$CB)
DimPlot(All, group.by = "inAlevin")
發(fā)現(xiàn)大部分的B細(xì)胞和T細(xì)胞沒(méi)有被鑒定出哩簿，而大部分的腫瘤細(xì)胞，漿細(xì)胞酝静，單核細(xì)胞被識(shí)別
這一點(diǎn)在alevin的issue (https://github.com/COMBINE-lab/salmon/issues/396) 中也提及：
"I think the CB frequency is most probably a bimodal distribution."
(可能是由于CB的雙相分布導(dǎo)致alevin對(duì)細(xì)胞CB的錯(cuò)誤估計(jì))
確實(shí)B和T細(xì)胞的文庫(kù)要偏小节榜，或許和細(xì)胞大小相關(guān)？
結(jié)論:
- alevin非潮鹬牵快宗苍，非常非常快薄榛。雖然實(shí)際分析還是用傳統(tǒng)比對(duì)流程好一些讳窟，但是快速的定量在多個(gè)數(shù)據(jù)集的探索，或者數(shù)據(jù)集質(zhì)控方面還是會(huì)有一些優(yōu)勢(shì)敞恋。
- alevin在比對(duì)的結(jié)果方面和Dropseq-Tools一致丽啡。
- alevin可能會(huì)錯(cuò)誤估計(jì)CB。
- 沒(méi)有看到alevin對(duì)傳統(tǒng)比對(duì)-定量流程的明顯優(yōu)勢(shì)硬猫，本來(lái)希望alevin在barcode和umi的錯(cuò)誤修復(fù)上會(huì)有優(yōu)勢(shì)补箍，目前看來(lái)優(yōu)勢(shì)不明顯。
- 目前還是傾向于選擇傳統(tǒng)的Dropseq-tools啸蜜，生成的bam文件同其他pipeline整合（如mutation calling）也是很方便的坑雅。
