UMI-tools count 使用

之前拿到了一組類(lèi)似Drop-seq的序列，雖然是雙端測(cè)序，但只有read2用作映射參考基因，read1只提供cell barcode+UMI用于后期去重。需要對(duì)序列匹配到同一個(gè)基因且擁有同樣cell barcode的唯一UMI進(jìn)行計(jì)數(shù)。于是就要用到UMI-tools的count命令。

Drop-seq的原理可以看文章 Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets

UMI-tools的使用說(shuō)明和教程可以在github上查看

1--安裝軟件

我比較習(xí)慣裝在conda里，還有其他的選擇可以從使用說(shuō)明里查看

$ 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 \ #設(shè)定read1里有UMI
--umi_len=8 \ #前8bp為UMI
--umi_prefix=UMI \ #命名
--trim_poly_x \ #去除read2 3'端的poly(A)尾
--thread=5 \ #線(xiàn)程數(shù)
-h ${cleandata}/${id}.fastp.html \ #結(jié)果報(bào)告html格式
-j ${cleandata}/${id}.fastp.json #結(jié)果報(bào)告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--映射到參考基因

因?yàn)閞ead1只用作提取UMI，這一步只用read2作為輸入

# hiast2
hisat2 -p 10 -x ${index} \ #構(gòu)建的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
輸出結(jié)果
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

由于需要根據(jù)匹配到的基因來(lái)去除UMI的重復(fù)，所以需要使用featureCounts將匹配到的基因標(biāo)記在序列中

# featureCounts
featureCounts \
-T 10 \ #線(xiàn)程
-g gene_id \ #從注釋文件中提取Meta-features信息用于read count，默認(rèn)是gene_id
-a ${gtf} \ #注釋的GTF文件
-o ${fcdir}/${id}.txt \ #輸出路徑和名字
-R BAM \ #輸出BAM文件(名字默認(rèn)為-輸入文件名.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
輸出結(jié)果

與上次輸出結(jié)果不同的地方是多了后三列，標(biāo)記了該序列是否匹配上基因（tag=XS）以及匹配到的具體基因編號(hào)（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

基本表達(dá)
umi_tools count [OPTIONS] --stdin=IN_BAM [--stdout=OUT_BAM]
部分參數(shù)說(shuō)明
參數(shù) 說(shuō)明
Barcode extraction options barcode提取選項(xiàng)
--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的標(biāo)簽
--cell-tag=CELL_TAG 包含cell barcode的標(biāo)簽
UMI grouping options UMI分組選項(xiàng)
--method=METHOD 用于UMI分組的方法。選擇：unique (讀取組共享完全相同的UMI); percentile; cluster; adjacency; directional [default] (確定連接的UMI的群集(基于漢明距離閾值)，并且umi A計(jì)數(shù)> =((2 * umi B計(jì)數(shù))-1)。每個(gè)網(wǎng)絡(luò)都是一個(gè)讀取組。)
single-cell RNA-Seq options 單細(xì)胞RNA-Seq選項(xiàng)
--per-gene 每個(gè)基因的組/重復(fù)/計(jì)數(shù)。必須與-gene-tag或--per-contig結(jié)合使用
--gene-tag=GENE_TAG 基因是由此bam標(biāo)簽定義的[default = none]
--assigned-status-tag=ASSIGNED_TAG Bam標(biāo)簽，描述是否將read指定給基因。默認(rèn)設(shè)置為與--gene-tag相同的標(biāo)簽
--per-cell 每個(gè)細(xì)胞的組/重復(fù)/計(jì)數(shù)
SAM/BAM options SAM/BAM文件選項(xiàng)
--mapping-quality=MAPPING_QUALITY 保留read的最低映射質(zhì)量[default= 0]
--unmapped-reads=UNMAPPED_READS 如何處理未映射的read。選擇：discard [default]; use; correct
--unpaired-reads=UNPAIRED_READS 如何處理未配對(duì)的read。選擇：discard; use [default]; correct
-i, --in-sam 輸入文件為sam格式 [default = False]
--paired 輸入BAM為雙端測(cè)序 [default = False]
-o, --out-sam 以sam格式輸出 [default = False]
--no-sort-output 不要對(duì)輸出進(jìn)行排序
input/output options 輸入/輸出選項(xiàng)
-I FILE, --stdin=FILE 輸入文件路徑 [default = stdin]
-S FILE, --stdout=FILE 輸出文件路徑 [default = stdout]

具體參數(shù)設(shè)置

umi_tools count \
--umi-separator=UMI_ \ #設(shè)置ID和UMI之間的分隔符為UMI_
--method=unique \ #UMI分組設(shè)置為unique，即只取唯一
--per-gene --gene-tag=XT \ #按基因計(jì)數(shù)，標(biāo)簽為XT
--assigned-status-tag=XS \ #設(shè)置匹配到基因的標(biāo)簽為XS
-I ${fcdir}/${id}.featureCounts.sorted.bam \ #輸入BAM文件路徑
-S ${umidir}/count/${id}.unique.counts.tsv.gz #輸出計(jì)數(shù)結(jié)果
輸出結(jié)果

為單條read的基因表達(dá)矩陣

gene                    count
ENSMUSG00000001138      1
ENSMUSG00000001143      6
ENSMUSG00000001674      1
ENSMUSG00000002881      1
ENSMUSG00000003134      11
ENSMUSG00000003135      4
ENSMUSG00000003458      1
ENSMUSG00000003464      2
ENSMUSG00000003721      3

如果要計(jì)算總和可以用下面的命令

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
輸出結(jié)果

多了末尾兩列，加上了UMI標(biāo)簽

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

結(jié)尾

由于我拿到的這組數(shù)據(jù)結(jié)果不是特別好，最后放棄了這組數(shù)據(jù)，也沒(méi)有進(jìn)行后續(xù)分析了，感覺(jué)UMI-tools還挺好用卻沒(méi)有教程，希望對(duì)大家有用！

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 今天感恩節(jié)哎，感謝一直在我身邊的親朋好友。感恩相遇！感恩不離不棄。 中午開(kāi)了第一次的黨會(huì)，身份的轉(zhuǎn)變要...
    余生動(dòng)聽(tīng)閱讀 10,798評(píng)論 0贊 11
  • 彩排完，天已黑
    劉凱書(shū)法閱讀 4,452評(píng)論 1贊 3
  • 沒(méi)事就多看看書(shū)，因?yàn)楦褂性?shī)書(shū)氣自華，讀書(shū)萬(wàn)卷始通神。沒(méi)事就多出去旅游，別因?yàn)闆](méi)錢(qián)而找借口，因?yàn)橹灰闶〕詢(xún)€用，來(lái)...
    向陽(yáng)之心閱讀 4,964評(píng)論 3贊 11
  • 表情是什么，我認(rèn)為表情就是表現(xiàn)出來(lái)的情緒。表情可以傳達(dá)很多信息。高興了當(dāng)然就笑了，難過(guò)就哭了。兩者是相互影響密不可...
    Persistenc_6aea閱讀 129,401評(píng)論 2贊 7

友情鏈接更多精彩內(nèi)容