古いDNAが断片化したサンプルからのシーケンスが増えている。しばしば数万年前のサンプルからも抽出される古代のサンプルのDNAは断片化が起きており、うまくDNAを抽出してもサイズが100-bpを超えることは滅多にない。短いDNAをペアードエンドでシーケンスすると、インサート全体がシーケンスされ、さらにアダプターまでシーケンスされることになる(図参照)。
ResearchGateより転載。
シーケンス全体がオーバーラップしているので、正しくマージできれば、単純に高いクオリテイのペア側の配列に更新するだけでシーケンス精度を高めることも可能になるが、元々のDNAの質が悪いので(ミスマッチとギャップが発生しうる)、マージは難しいと推測される。
leeHomは5'側と3'側のアダプターを除去し、ベイジアン最大事後確率アプローチを用いて元のDNA配列を再構成する方法論となる。シミュレーションと、古代のサンプルとして有名なネアンデルタール人のシーケンスデータを使いテストされており、他の方法論より精度が高いという結果が出ている。シングルエンド、ペアードエンドのシーケンスデータに対応している。
公式サイト
https://bioinf.eva.mpg.de/leehom/
インストール
https://github.com/grenaud/leeHom
git clone --recursive https://github.com/grenaud/leeHom.git
cd leeHom/
make
src/leeHom #動作確認
$ src/leeHom
Usage:
src/leeHom [options] BAMfile
This program takes an unaligned BAM where mates are consecutive
or fastq files and trims and merges reads
You can specify a unaligned bam file or one or two fastq :
-fq1 First fastq
-fq2 Second fastq file (for paired-end)
-fqo Output fastq prefix
-o , --outfile Output (BAM format).
-u Produce uncompressed bam (good for pipe)
--aligned Allow reads to be aligned (default false)
-v , --verbose Turn all messages on (default false)
--log [log file] Print a tally of merged reads to this log file (default only to stderr)
--phred64 Use PHRED 64 as the offset for QC scores (default : PHRED33)
Paired End merging/Single Read trimming options
You can specify either:
--ancientdna ancient DNA (default false)
this allows for partial overlap
or if you know your size length distribution:
--loc Location for lognormal dist. (default none)
--scale Scale for lognormal dist. (default none)
--keepOrig Write original reads if they are trimmed or merged (default false)
Such reads will be marked as PCR duplicates
-f , --adapterFirstRead Adapter that is observed after the forward read (def. Multiplex: AGATCGGAAGAGCACACGTCTGAACTCCAG)
-s , --adapterSecondRead Adapter that is observed after the reverse read (def. Multiplex: AGATCGGAAGAGCGTCGTGTAGGGAAAGAG)
-c , --FirstReadChimeraFilter If the forward read looks like this sequence, the cluster is filtered out.
Provide several sequences separated by comma (def. Multiplex: ACACTCTTTCCCTACACGTCTGAACTCCAG)
-k , --key Key sequence with which each sequence starts. Comma separate for forward and reverse reads. (default '')
-i , --allowMissing Allow one base in one key to be missing or wrong. (default false)
--trimCutoff Lowest number of adapter bases to be observed for single Read trimming (default 1)
パスを通しておく。
マルチコア版
> src/leeHomMulti #-tでthread数を決める
$ src/leeHomMulti
Usage:
leeHomMulti [options] BAMfile
This program takes an unaligned BAM where mates are consecutive
or fastq files and trims and merges reads
You can specify a unaligned bam file or one or two fastq :
-fq1 First fastq
-fq2 Second fastq file (for paired-end)
-fqo Output fastq prefix
-o , --outfile Output (BAM format).
-u Produce uncompressed bam (good for pipe)
--aligned Allow reads to be aligned (default false)
-v , --verbose Turn all messages on (default false)
--log [log file] Print a tally of merged reads to this log file (default only to stderr)
--phred64 Use PHRED 64 as the offset for QC scores (default : PHRED33)
-t [threads] Use multiple cores (default : 1)
Paired End merging/Single Read trimming options
You can specify either:
--ancientdna ancient DNA (default false)
this allows for partial overlap
or if you know your size length distribution:
--loc Location for lognormal dist. (default none)
--scale Scale for lognormal dist. (default none)
--keepOrig Write original reads if they are trimmed or merged (default false)
Such reads will be marked as PCR duplicates
-f , --adapterFirstRead Adapter that is observed after the forward read (def. Multiplex: AGATCGGAAGAGCACACGTCTGAACTCCAG)
-s , --adapterSecondRead Adapter that is observed after the reverse read (def. Multiplex: AGATCGGAAGAGCGTCGTGTAGGGAAAGAG)
-c , --FirstReadChimeraFilter If the forward read looks like this sequence, the cluster is filtered out.
Provide several sequences separated by comma (def. Multiplex: ACACTCTTTCCCTACACGTCTGAACTCCAG)
-k , --key Key sequence with which each sequence starts. Comma separate for forward and reverse reads. (default '')
-i , --allowMissing Allow one base in one key to be missing or wrong. (default false)
--trimCutoff Lowest number of adapter bases to be observed for single Read trimming (default 1)
ラン
テストラン。
fasrqファイルからアダプターを除く。
leeHom -f AGATCGGAAGAGCACACGTCTGAACTCCAG -s GGAAGAGCGTCGTGTAGGGAAAGAGTGTAG --ancientdna -fq1 testData/rawAncientDNA.f1.gz -fq2 testData/rawAncientDNA.f2.gz -fqo testData/outfq
- -f Adapter that is observed after the forward read (def. Multiplex: AGATCGGAAGAGCACACGTCTGAACTCCAG)
- -s Adapter that is observed after the reverse read (def. Multiplex: AGATCGGAAGAGCGTCGTGTAGGGAAAGAG)
- --ancientdna ancient DNA (default false) this allows for partial overlap
- -fq1 First fastq-fq1 First fastq
- -fq2 Second fastq file (for paired-end)
- -fqo Output fastq prefix
ランが終わると下の解析logが表示される。
Total 50000; Merged (trimming) 40540; Merged (overlap) 7376; Kept PE/SR 1955; Trimmed SR 0; Adapter dimers/chimeras 129; Failed Key 0
47,916のペアードエンドがマージされ(=40540+7376)、1955は重複が見つからずマージされなかったことを意味している。またインサートなしにアダプターが2つタンデムに結合したキメラは129見つかった。
出力はマージされたfastq.gzと、マージされなかったr1.fastq.gz、r2.fastq.gz、そして条件が満たされなかったfailのfastq.gzが出力される。
bamファイルからアダプターを除く。
leeHom -f AGATCGGAAGAGCACACGTCTGAACTCCAG -s GGAAGAGCGTCGTGTAGGGAAAGAGTGTAG --ancientdna -o testData/reconsAncientDNA.bam testData/rawAncientDNA.bam
- -o , --outfile Output (BAM format).
引用
leeHom: adaptor trimming and merging for Illumina sequencing reads
Gabriel Renaud, Udo Stenzel, and Janet Kelso.
Nucleic Acids Res. 2014 Oct 13; 42(18): e141.