一、需求背景
10X V5轉(zhuǎn)錄組vdj分析中，因?yàn)橛?strong>一段特殊基因序列比對不上參考基因組，所以bam文件沒有這個(gè)barcode信息。需要把R2的這段序列，根據(jù)序列ID從R1中找出R2序列與R1 barcode和UMI的對應(yīng)關(guān)系，然后做UMI去重，統(tǒng)計(jì)數(shù)量，然后再映射到細(xì)胞聚類圖中去看看這個(gè)序列的UMI表達(dá)量
二、思路分析
根據(jù)文庫結(jié)構(gòu)
圖片來源：https://blog.csdn.net/herokoking/article/details/103629141
將目標(biāo)序列作為比對參考模板，用R2文件進(jìn)行blast比對，從比對結(jié)果中篩選出從3'到5’方向的，因?yàn)?0X V5轉(zhuǎn)錄組的R2的測序方向就是3’到5’，然后再篩選出比對上序列長度大于等于85nt的子集，獲得比對上fastq R2測序文件的reads的ID信息， 再通過seqtk工具查詢在 fsatq R1文件根據(jù)ID提取出相應(yīng)的序列信息。
從查詢出來的序列中提取序列前1到16nt獲取barcode，17到28nt獲得UMI信息，經(jīng)過barcode和UMI的序列去重、統(tǒng)計(jì)、分群的標(biāo)簽映射
流程圖：
序列模板---->R2 blast比對----->過濾&提取ID---->seqtk提取R1序列--->barcode和UMI的序列去重-->數(shù)據(jù)統(tǒng)計(jì)--->分群的標(biāo)簽映射---->UMI表達(dá)量可視化分析
三、代碼
3.1 blast比對
#fastq轉(zhuǎn)換為fasta
nohup less -S ../sample-*_L3_2.fq.gz |awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' > ./sample-*_L3_2.fasta &
# 數(shù)據(jù)庫構(gòu)建
makeblastdb \
-dbtype nucl \
-in artificial.fasta \
-input_type fasta \
-parse_seqids \
-out artificial.blastdb
# 比對
##這里的*號(hào)表示4條R2序列的編號(hào)省略
blastn \
-query sample-*_L3_2.fasta \
-db artificial.blastdb \
-out sample-*_results2.xls \
-outfmt 6 \
-num_threads 8
3.2 過濾比對上的reads2序列
cat ../sample-1_results2.xls |sort -nr -k4|awk '$4>=85 && $9>$10{print $0}' > map_sample-1.xls
cat ../sample-2_results2.xls |sort -nr -k4|awk '$4>=85 && $9>$10{print $0}' > map_sample-2.xls
cat ../sample-3_results2.xls |sort -nr -k4|awk '$4>=85 && $9>$10{print $0}' > map_sample-3.xls
cat ../sample-4_results2.xls |sort -nr -k4|awk '$4>=85 && $9>$10{print $0}' > map_sample-4.xls
3.3 R2 序列id獲取
ls *.xls |cat |xargs awk '{print $1}' >query_ids2
3.3 合并R1序列方便后續(xù)查找操作
ls sample*_1.fq.gz |xargs less >sample_R1.fastq
3.4 使用seqtk工具根據(jù)R2 ID提取R1序列子集
seqtk安裝、部署、使用方法參考：
# 輸出fastq格式
$seqtk subseq sample_R1.fastq query_ids2 >seqtk_result.tsv
# -t 參數(shù)：輸出一種\t分隔格式，更方便后面的提取barcode操作
$seqtk subseq -t sample_R1.fastq query_ids2 >seqtk_result2.tsv
3.5 barcode去重，用于生成有無表達(dá)特定序列的分組信息標(biāo)簽
cat seqtk_result2.tsv |cut -f3 |awk '{print substr($1, 1, 28)}' |sort |uniq |awk '{print substr($1,1,16)}' |sort |uniq >uniq_barcode.tsv
3.6 與關(guān)注細(xì)胞細(xì)分亞群的barcode進(jìn)行取交集操作，過濾掉無關(guān)的序列
R:
cellType<-readRDS('cellType.rds')
write.table(gsub('-.*', '', rownames(cellType@meta.data)), file = 'cellType_barcode.tsv', sep = '\t', quote = F, col.names = F, row.names = F)
Linux:
# 取交集
sort cellType_barcode.tsv uniq_barcode.tsv |uniq -d >cellType_intersect.tsv
3.7 根據(jù)barcode和UMI一起去重
cat seqtk_result2.tsv |cut -f3 |awk '{print substr($1, 1, 28)}' |sort |uniq >uniq_barcode_umi.tsv
3.8 統(tǒng)計(jì)每個(gè)細(xì)胞barcode相應(yīng)的關(guān)注序列的UMI
cat uniq_barcode_umi.tsv |awk '{hash[substr($1, 1, 16)]+=1}END{for (i in hash){printf("%s\t%d\n", i, hash[i])}}' >count_umi.tsv
3.9 將分組信息和UMI信息映射到關(guān)注的分群聚類圖中（數(shù)據(jù)可視化）
R:
library(plyr)
library(Seurat)
library(ggplot2)
#read barcode group label
interBarcodes<-read.table('cellType_intersect.tsv', sep = '\t', header = F)
#read UMI label
umi_count<-read.table('count_umi.tsv', sep='\t', header = F, col.names = c('barcode', 'umi'))
#mapping group label
cellType$artificial<-ifelse(gsub('-.*', '', rownames(cellType@meta.data))%in%interBarcodes$V1,
'artificial_gene', 'no_artificial_gene')
#mapping UMI label
cellType$artificial_gene<-mapvalues(gsub('-.*', '', rownames(cellType@meta.data)), as.character(umi_count$barcode), as.character(umi_count$umi))
cellType$artificial_gene<-as.numeric(cellType$artificial_gene)
#change the no expression UMI label from NA to 0
cellType$artificial_gene<-ifelse(is.na(cellType$artificial_gene), 0, cellType$artificial_gene)
head(cellType@meta.data)
# Visualize
DimPlot(cellType, group.by = 'artificial', reduction = 'tsne')
FeaturePlot(cellType, features = 'artificial_gene', reduction = 'tsne', label=T, cols = c('grey','red'), pt.size = 1)
