长读长序列快速比对工具Minimap2实战之cs tag

作为长读长测序快速比对的经典工具，Minimap2目前已更新至2.24版本。其基本用法已被大家熟知：

git clone https://github.com/lh3/minimap2
cd minimap2 && make
# long sequences against a reference genome
./minimap2 -a test/MT-human.fa test/MT-orang.fa > test.sam
# create an index first and then map
./minimap2 -x map-ont -d MT-human-ont.mmi test/MT-human.fa
./minimap2 -a MT-human-ont.mmi test/MT-orang.fa > test.sam
# use presets (no test data)
./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio CLR genomic reads
./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
./minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or later)
./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio HiFi/CCS genomic reads (v2.18 or earlier)
./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # Final PacBio Iso-seq or traditional cDNA
./minimap2 -ax splice --junc-bed anno.bed12 ref.fa query.fa > aln.sam  # prioritize on annotated junctions
./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
# man page for detailed command line options
man ./minimap2.1

但需要注意的是，由于Minimap2针对的是长读长的三代测序数据，其产生的比对结果可能会比较复杂，体现在CIGAR上就是操作符可能超过了65535个，这显然超过了bam文件能够接受的上限，比如使用Picard或Samtools将sam文件转为bam文件时，就会报错。这是就需要在调用Minimap2时加上"-L"参数，该参数将CIGAR转移至CG标签上，这样进行下游操作时就不会报错。

再比如，通过cs标签可以记录query序列的哪些碱基出现了SNP和INDEL，其通配格式如下：

/(:[0-9]+|\*[a-z][a-z]|[=\+\-][A-Za-z]+)+/

cs标签是CIGAR的细化，打头的字符代表相应的CIGAR操作，紧随其后的就是相应的序列，cs标签需要"--cs"参数来打开，默认为"--cs=short"。如下比对结果其cs标签为:6-ata:10+gtc:4*at:3，:[0-9]+ 代表完全匹配, -ata代表缺失, +gtc为插入，*at代表参考碱基a被query碱基t替代。当使用"--cs=long"参数时，cs标签会将完全匹配的碱基也显示出来，如=CGATCG-ata=AATAGAGTAG+gtc=GAAT*at=GCA。

CGATCGATAAATAGAGTAG---GAATAGCA
||||||   ||||||||||   |||| |||
CGATCG---AATAGAGTAGGTCGAATtGCA

为什么会提到cs标签呢？因为某些软件在解析bam或sam文件时，会读取cs标签里的内容，否则会报错，比如cs标签未找到。大家在使用过程中，可以灵活把握，下游分析软件报错时，检查一下是否是CIGAR或者是cs标签的问题。