之前拿到了一组类似Drop-seq的序列，虽然是双端测序，但只有read2用作映射参考基因，read1只提供cell barcode+UMI用于后期去重。需要对序列匹配到同一个基因且拥有同样cell barcode的唯一UMI进行计数。于是就要用到UMI-tools的count命令。
Drop-seq的原理可以看文章 Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets
1--安装软件
我比较习惯装在conda里，还有其他的选择可以从使用说明里查看
$ conda install -c bioconda -c conda-forge umi_tools
2--在序列head里加上UMI信息
首先查看原始序列，read1的前8bp是UMI，read2末尾有一段poly(A)尾需要去除
原始序列
## read1
@A00159:637:HC555DSXY:2:1443:4679:29387 1:N:0:CTAGGAAT+TATAGCCT
CGCTCGCCTCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTAAAATCAATTTAAGTTTGTGATGAAGATTTGGATATTTAGTGAGGAAAAAGAGAAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:,:,,,,,,,,,,,,,,,,,:,,,:,,,,,,,,,,,,,,,,,,,,,:,,,,,,,,,,,,,,,,,,,
## read2
@A00159:637:HC555DSXY:2:1443:4679:29387 2:N:0:CTAGGAAT+TATAGCCT
GGACTGGATCCCGTCCTCTAACGTCCGTCTCTGGCCGAGTCTAACACTGTACAACTGTCTCTGACCATTAAATGCTGTTGTACCGTGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF:FFFF,FFF:F:F::F::FFF,:
使用fastp添加UMI
fastp \
-i ${fqdir}/${id}_1.clean.fq.gz \ #原始read1
-I ${fqdir}/${id}_2.clean.fq.gz \ #原始read2
-o ${fpdir}/${id}_1.fastp.fq.gz \ #输出read1
-O ${fpdir}/${id}_2.fastp.fq.gz \ #输出read2
-Q \ #disable quality filtering
-U \ #提取UMI
--umi_loc=read1 \ #设定read1里有UMI
--umi_len=8 \ #前8bp为UMI
--umi_prefix=UMI \ #命名
--trim_poly_x \ #去除read2 3'端的poly(A)尾
--thread=5 \ #线程数
-h ${cleandata}/${id}.fastp.html \ #结果报告html格式
-j ${cleandata}/${id}.fastp.json #结果报告json格式
处理后的序列
# read1
@A00159:637:HC555DSXY:2:1443:4679:29387:UMI_CGCTCGCC 1:N:0:CTAGGAAT+TATAGCCT
TCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTAAAATCAATTTAAGTTTGTGATGAAGATTTGGATATTTAGTGAGGAAAAAGAGAAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:,:,,,,,,,,,,,,,,,,,:,,,:,,,,,,,,,,,,,,,,,,,,,:,,,,,,,,,,,,,,,,,,,
# read2
@A00159:637:HC555DSXY:2:1443:4679:29387:UMI_CGCTCGCC 2:N:0:CTAGGAAT+TATAGCCT
GGACTGGATCCCGTCCTCTAACGTCCGTCTCTGGCCGAGTCTAACACTGTACAACTGTCTCTGACCATTAAATGCTGTTGTACCGTGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
3--映射到参考基因
因为read1只用作提取UMI，这一步只用read2作为输入
# hiast2
hisat2 -p 10 -x ${index} \ #构建的index
-U ${fpdir}/read1/${id}_2.fastp.fq.gz \ #只输入read2
-S ${hisdir}/read1/${id}.hisat2.sam #输出SAM文件
# samtools sort & index
samtools view -bS ${hisdir}/${id}.hisat2.sam | \
samtools sort -@ 10 -o ${samdir}/${id}.hisat2.sorted.bam && \
samtools index ${samdir}/${id}.hisat2.sorted.bam ${samdir}/${id}.hisat2.sorted.bam.bai
输出结果
A00159:637:HC555DSXY:2:1443:4679:29387:UMI_CGCTCGCC 0 chr10 127281719 60 88M * 0 0 GGACTGGATCCCGTCCTCTAACGTCCGTCTCTGGCCGAGTCTAACACTGTACAACTGTCTCTGACCATTAAATGCTGTTGTACCGTGG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:-5 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:25T62 YT:Z:UU NH:i:1
4--featureCounts
由于需要根据匹配到的基因来去除UMI的重复，所以需要使用featureCounts将匹配到的基因标记在序列中
# featureCounts
featureCounts \
-T 10 \ #线程
-g gene_id \ #从注释文件中提取Meta-features信息用于read count，默认是gene_id
-a ${gtf} \ #注释的GTF文件
-o ${fcdir}/${id}.txt \ #输出路径和名字
-R BAM \ #输出BAM文件(名字默认为-输入文件名.featureCounts.bam)
${samdir}/${id}.hisat2.sorted.bam #输入BAM文件
# samtools sort & index
samtools sort ${fcdir}/${id}.hisat2.sorted.bam.featureCounts.bam -o ${fcdir}/${id}.featureCounts.sorted.bam && \
samtools index ${fcdir}/${id}.featureCounts.sorted.bam
输出结果
与上次输出结果不同的地方是多了后三列，标记了该序列是否匹配上基因（tag=XS）以及匹配到的具体基因编号（tag=XT）
A00159:637:HC555DSXY:2:1443:4679:29387:UMI_CGCTCGCC 0 chr10 127281719 60 88M * 0 0 GGACTGGATCCCGTCCTCTAACGTCCGTCTCTGGCCGAGTCTAACACTGTACAACTGTCTCTGACCATTAAATGCTGTTGTACCGTGG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:-5 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:25T62 YT:Z:UU NH:i:1 XS:Z:Assigned XN:i:1 XT:Z:ENSMUSG00000025410
5--UMI-tools count
基本表达
umi_tools count [OPTIONS] --stdin=IN_BAM [--stdout=OUT_BAM]
部分参数说明
| 参数 | 说明 |
|---|---|
| Barcode extraction options | barcode提取选项 |
| --extract-umi-method=GET_UMI_METHOD | cell barcode和UMI是如何编码？选择：read_id [default]; tag; umis |
| --umi-separator=UMI_SEP | 读取ID和UMI之间的分隔符 |
| --umi-tag=UMI_TAG | 包含UMI的标签 |
| --cell-tag=CELL_TAG | 包含cell barcode的标签 |
| UMI grouping options | UMI分组选项 |
| --method=METHOD | 用于UMI分组的方法。选择：unique (读取组共享完全相同的UMI); percentile; cluster; adjacency; directional [default] (确定连接的UMI的群集(基于汉明距离阈值)，并且umi A计数> =((2 * umi B计数)-1)。每个网络都是一个读取组。) |
| single-cell RNA-Seq options | 单细胞RNA-Seq选项 |
| --per-gene | 每个基因的组/重复/计数。必须与-gene-tag或--per-contig结合使用 |
| --gene-tag=GENE_TAG | 基因是由此bam标签定义的[default = none] |
| --assigned-status-tag=ASSIGNED_TAG | Bam标签，描述是否将read指定给基因。默认设置为与--gene-tag相同的标签 |
| --per-cell | 每个细胞的组/重复/计数 |
| SAM/BAM options | SAM/BAM文件选项 |
| --mapping-quality=MAPPING_QUALITY | 保留read的最低映射质量[default= 0] |
| --unmapped-reads=UNMAPPED_READS | 如何处理未映射的read。选择：discard [default]; use; correct |
| --unpaired-reads=UNPAIRED_READS | 如何处理未配对的read。选择：discard; use [default]; correct |
| -i, --in-sam | 输入文件为sam格式 [default = False] |
| --paired | 输入BAM为双端测序 [default = False] |
| -o, --out-sam | 以sam格式输出 [default = False] |
| --no-sort-output | 不要对输出进行排序 |
| input/output options | 输入/输出选项 |
| -I FILE, --stdin=FILE | 输入文件路径 [default = stdin] |
| -S FILE, --stdout=FILE | 输出文件路径 [default = stdout] |
具体参数设置
umi_tools count \
--umi-separator=UMI_ \ #设置ID和UMI之间的分隔符为UMI_
--method=unique \ #UMI分组设置为unique，即只取唯一
--per-gene --gene-tag=XT \ #按基因计数，标签为XT
--assigned-status-tag=XS \ #设置匹配到基因的标签为XS
-I ${fcdir}/${id}.featureCounts.sorted.bam \ #输入BAM文件路径
-S ${umidir}/count/${id}.unique.counts.tsv.gz #输出计数结果
输出结果
为单条read的基因表达矩阵
gene count
ENSMUSG00000001138 1
ENSMUSG00000001143 6
ENSMUSG00000001674 1
ENSMUSG00000002881 1
ENSMUSG00000003134 11
ENSMUSG00000003135 4
ENSMUSG00000003458 1
ENSMUSG00000003464 2
ENSMUSG00000003721 3
如果要计算总和可以用下面的命令
zcat ${id}.unique.counts.tsv.gz | perl -alne '$sum += $F[1]; END {print $sum}'
如果需要输出的BAM文件，可以用UMI-tools中的group命令
umi_tools group \
--umi-separator=UMI_ \
--method=unique \
--per-gene --gene-tag=XT \
-I ${fcdir}/${id}.featureCounts.sorted.bam \
--group-out=${umidir}/${id}.grouped.tsv \
--log=${umidir}/${id}.group.log \
--output-bam | \
samtools view -o ${umidir}/${id}.umi.bam
输出结果
多了末尾两列，加上了UMI标签
A00159:637:HC555DSXY:2:1443:4679:29387:UMI_CGCTCGCC 0 chr10 127281719 60 88M * 0 0 GGACTGGATCCCGTCCTCTAACGTCCGTCTCTGGCCGAGTCTAACACTGTACAACTGTCTCTGACCATTAAATGCTGTTGTACCGTGG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:-5 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:25T62 YT:Z:UU NH:i:1 XS:Z:Assigned XN:i:1 XT:Z:ENSMUSG00000025410 UG:i:4183 BX:Z:CGCTCGCC
结尾
由于我拿到的这组数据结果不是特别好，最后放弃了这组数据，也没有进行后续分析了，感觉UMI-tools还挺好用却没有教程，希望对大家有用！
