Githubより
ほとんどのショートリードアライナーは、参照ゲノムに対するリードのローカルアライメントを行う。例として、bwa mem、minimap2、bowtie2などがある（--end-to-endモードの場合を除く）。つまり、リードの両端がベストアライメントに含まれていない可能性がある。これは次のような場合に起こる。
- アダプター配列（リファレンスにはない）
- クオリティの低い塩基（ミスマッチがあるとアライメントのスコアが悪くなる）
- リファレンスと比較したサンプルの構造上の差異
- コンティグの最初と最後のリード（環状ゲノムを含む）
リードアライナーはSAMファイルを出力する。このフォーマットの6列目には、CIGAR文字列が格納されており、リードのどの部分がアラインメントされ、どの部分がアラインメントされなかったかが記述されている。アラインメントされなかったリードの端は、「ソフト」または「ハード」にクリップされ、CIGAR文字列の両端にSとHで示される。両方のタイプが存在する可能性もあるが、一般的ではない。ソフトクリップとハードクリップは生物学的には何の意味もなく、SAMファイルに完全な読み取り配列があるかどうかを指しているだけである。
100bpのアラインドリードの例。
- 100M - 完全にアラインメントされたリード
- 30S70M - 左端の30bpをソフトクリップした後、70bpをアラインしたもの。
- 95M5H - 95bpのアラインメントと、それに続くアラインメントされていない5bp（ハードクリップ）。
- 25S40M35S - リードの真ん中の40bpだけがアラインされている。
動機
下流の問題を避けるために、これらのアラインメントを除去したい場合がある。特にsamclipは、偽陽性の原因となる疑わしいクリップされたアラインメント部分を除去することで、バリアントコーリングを改善するように設計されている。ただし、コンティグの端に位置する場合には、そのアラインメントは維持される。
インストール
#conda (link)
mamba install -c bioconda -c conda-forge samclip
#brew
brew install brewsci/bio/samclip
> ./samclip --help
$ samclip --help
SYNOPSIS
Filter SAM file for soft & hard clipped alignments
AUTHOR
Torsten Seemann (@torstenseemann)
USAGE
% samclip --ref ref.fa < in.sam > out.sam
% minimap2 ref.fa R1.fq R2.fq | samclip --ref ref.fa | samtools sort > out.bam
OPTIONS
--help This help
--version Print version and exit
--ref FASTA Reference genome - needs FASTA.fai index
--max NUM Maximum clip length to allow (default=5)
--invert Output rejected SAM lines and ignore good ones
--debug Print verbose debug info to stderr
--progress N Print progress every NUM records (default=100000,none=0)
HOMEPAGE
https://github.com/tseemann/samclip
実行方法
samファイルを指定する。
#sam
samclip --ref ref.fa < in.sam > out.sam
#bam
samtools view -h in.bam | samclip --ref ref.fa | samtools sort > out.bam
#fastq
bwa mem ref.fa R1.fq R2.fq | samclip --ref ref.fa | samtools sort > out.bam
クリップされたリードを出力する。
#sam
samclip --invert --ref ref.fa < in.sam > clipped.sam
引用
https://github.com/tseemann/samclip
