macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

コンタミやダメージを考慮してAncient DNAのシーケンシングリードをシミュレートする gargammel

 

 Ancient DNA(aDNA)とも呼ばれるsubfossilsから回収されたDNAは、populationの歴史を再構築するためにますます使用されている(Leonardi et al、2016)。しかし、下流の推論に影響を与える可能性があるいくつかの要因があるため、aDNAデータの分析は依然として困難である。第一に、DNAは時間とともに分解する傾向があり、かなりのヌクレオチドのmisincorporationsを示す限られたサイズ(30〜80 bp)のフラグメントをもたらす(Briggs et al、2007)。第二に、生物死後に環境微生物がコロニーを形成する傾向がある(Green et al、2009 pubmed)。その結果、内因性DNA比率が極端に減少することがあり、ショットガンシークエンシングアプローチが非経済的になる。第三に、そのような外因性配列は、リードアラインメント中に正しく同定されないと、Ancient DNAの再構築に影響を及ぼし得る。ヒト族から回収されたDNAの場合、発掘中および実験室中を含む任意の段階で導入され得る現在のヒトDNA混入は、単一のサンプル内に無関係な集団の歴史を混在させるので特に問題がある。

 Ancient DNA研究者は、人口パラメータの推論を目的とした要約統計量の頑健性をテストするためにシミュレーションを頻繁に使用する。いくつかのパッケージはシーケンシングに起因するプラットフォーム特有のエラーをシミュレートしているが、損傷、断片化、ヒトおよび微生物汚染などのそれらの最も顕著な特徴を含むaDNA配列データセットを適切にシミュレートするパッケージは現在利用できない。

 ここでは、微生物画分、内因性画分、そして今日の人間の汚染を表す一連のゲノム参照からのaDNA配列データセットをシミュレートするパッケージであるgargammelを紹介する。このパッケージは、死後のDNA損傷や塩基の誤取り込みなど、DNA配列の最も一般的な機能をシミュレートすることができる。さらに、ライブラリーの調製に使用される分子ツールからのバイアス、GCに富むフラグメントからのバイアス、およびシーケンシングプラットフォームによってもたらされるエラーをシミュレートする。

 

 

f:id:kazumaxneo:20190315104433j:plain

Gargammel flowchart.  論文より転載。

 

HP

gargammel by grenaud

f:id:kazumaxneo:20190315110944j:plain

 


インストール

依存

  • git
  • C++ compiler supporting C++11
  • cmake, you can install on Ubuntu by typing: sudo apt install cmake
  • zlib
  • lib gsl, you can install on Ubuntu by typing: sudo apt-get install libgsl0-dev

If you plan on using ms2chromosomes.py to simulate chromosomes based on ms, you also need:

本体 Github

git clone https://github.com/grenaud/gargammel.git
cd gargammel/
make

perl gargammel.pl

# perl gargammel.pl 

 

 

 

 This script is a wrapper to run the different programs to create

 a set of Illumina reads for ancient DNA from a set of fasta files

 representing the endogenous, the contaminant from the same species

 and the bacterial contamination.

 

 

 

 usage: gargammel.pl <options> [directory with fasta directories] 

 

 This directory should contain 3 directories:

  bact/ The bacterial contamination

  cont/ The contamination from the same species

  endo/ The endogenous material

 

 Options:

--comp [B,C,E] Composition of the final set in fraction 

the 3 numbers represent the bacterial, contaminant and endogenous

ex: --comp 0.6,0.02,0.38 will result

in 60% bacterial contamination while the rest will be from the same

species 5% will be contamination and 95% will be endogenous

Default: --comp 0,0,1

--mock Do nothing, just print the commands that will be run

-o Output prefix (default: [input dir]/simadna)

 Either specify:

-n [number] Generate [number] fragments (default: 1000)

-c [coverage] Endogenous coverage

 

 Fragment selection

 ===================

  --misince [file] Base misincorporation for the endogenous fragments (default none)

  --misincc [file] Base misincorporation for the contaminant fragments (default none)

  --misincb [file] Base misincorporation for the bacterial fragments (default none)

  --distmis [dist] Distance to consider for base misincorporation (default 1)

  this file is obtained from mapDamage

 

  Fragment size distribution: specify either one of the 3 possible options:

-l [length] Generate fragments of fixed length  (default: 35)

-s [file] Open file with size distribution (one fragment length per line)

-f [file] Open file with size frequency in the following format:

    length[TAB]freq ex:

    40 0.0525

    41 0.0491

    ...

 

Length options:

--loc [file] Location for lognormal distribution (default none)

--scale [file] Scale for lognormal distribution    (default none)

 

  Fragment size limit:

--minsize [length] Minimum fragments length (default: 0)

                --maxsize [length] Maximum fragments length (default: 1000)

 

  Fragment methylation:

--methyl Allow for lowercase C and G to denote

                methylated cytosines on the + and - strand respectively

                (default: not used)

 

 Deamination

 ===================

 To add deamination to the bacterial and endogenous material,you can specify

 either one of these options:

-mapdamage [mis.txt] [protocol] Read the miscorporation file [mis.txt]

                                  produced by mapDamage

                      [protocol] can be either "single" or "double" (without quotes)

                      Single strand will have C->T damage on both ends

                      Double strand will have and C->T at the 5' end and G->A damage at the 3' end

-matfile [matrix file prefix] Read the matrix file of substitutions

                                Provide the prefix only, both files must end with

                              5.dat and 3.dat

-damage [v,l,d,s]   For the Briggs et al. 2007 model

                  The parameters must be comma-separated e.g.: -damage 0.03,0.4,0.01,0.3

                  v: nick frequency

                  l: length of overhanging ends (geometric parameter)

                  d: prob. of deamination of Cs in double-stranded parts

                  s: prob. of deamination of Cs in single-stranded parts

 

 Alternatively, you can specify these options independently for the endogenous (e), bacterial (b)

 and present-day human contaminant (c) using the following options:

-mapdamagee [mis.txt] [protocol] Endogenous mapDamage misincorporation file

-matfilee [matrix file prefix] Endogenous matrix file of substitutions

-damagee [v,l,d,s]   Endogenous Briggs parameters

-mapdamageb [mis.txt] [protocol] Bacterial mapDamage misincorporation file

-matfileb [matrix file prefix] Bacterial matrix file of substitutions

-damageb [v,l,d,s]   Bacterial Briggs parameters

-mapdamagec [mis.txt] [protocol] Human contaminant mapDamage misincorporation file

-matfilec [matrix file prefix] Human contaminant matrix file of substitutions

-damagecd [v,l,d,s]   Human contaminant Briggs parameters

 

 please note that if you do specify deamination for one source but not for another, no deamination will be added

 

 If using --methyl, you can also specify different matrix file for methylated

-matfilenonmeth [matrix file prefix] Read the matrix file of substitutions for non-methylated Cs

                            Provide the prefix only, both files must end with

                            5.dat and 3.dat

   -matfilemeth [matrix file prefix] Read the matrix file of substitutions for methylated Cs

                        Provide the prefix only, both files must end with

                        5.dat and 3.data

 Adapter and sequencing

 ===================

-fa [seq] Adapter that is observed after the forward read (Default: AGATCGGAAG...)

-sa [seq] Adapter that is observed after the reverse read (Default: AGATCGGAAG...)

-rl [length] Desired read length  (Default: 75)

-se                             use single-end sequencing (Default: paired-end)

 

                                        The following options change the sequencing error rate, please note that positive factor

                                        will decrease the rate of such errors and a negative one will increase it.

        -qs     [factor]                Increase error rate for forward reads by a factor of 1/(10^([factor]/10)) (Default: 0)

        -qs2    [factor]                Increase error rate for reverse reads by a factor of 1/(10^([factor]/10)) (Default: 0)

 

 

-ss     [system]                Illumina platfrom to use, the parentheses indicate the max. read length

                                use the shorthand in the left column:

                                                                        (single-end, paired-end)

  GA2  - GenomeAnalyzer II (  50bp,  75bp)

  HS20 - HiSeq 2000        ( 100bp, 100bp)

  HS25 - HiSeq 2500        ( 125bp, 150bp) (Default)

  HSXt - HiSeqX TruSeq     ( 150bp, 150bp)

  MSv1 - MiSeq v1          ( 250bp, 250bp)

  MSv3 - MiSeq v3          ( 250bp, 250bp)

 

 

 

 

依存ツールも含めてビルドされる。しかしsamtoolsがビルドされなかったので、手っ取り早くapt installで古いバージョンを入れた。

 

dockerイメージも上げておきます。

$ docker run -itv $PWD:/data/ kazumax/gargammel
> source ~/.profile
> gargammel.pl -h

 

テストラン

1、リファレンスを作る。

cd gargammel/ && mkdir data && cd data/

#実行
python ../ms2chromosomes.py -s 0.2 -f . -n 1000

#cleanup(作業ファイルが)
rm -rfv simul_* seedms

ls -alh *

# ls -alh *

-rw-r--r-- 1 root root 9.6M Mar 15 11:01 ref.fa

 

bact:

total 132K

drwxr-xr-x 2 root root 4.0K Mar 15 11:01 .

drwxr-xr-x 5 root root 124K Mar 15 11:02 ..

 

cont:

total 20M

drwxr-xr-x 2 root root 4.0K Mar 15 11:01 .

drwxr-xr-x 5 root root 124K Mar 15 11:02 ..

-rw-r--r-- 1 root root 9.6M Mar 15 11:01 cont.1.fa

-rw-r--r-- 1 root root 9.6M Mar 15 11:01 cont.2.fa

 

endo:

total 20M

drwxr-xr-x 2 root root 4.0K Mar 15 11:01 .

drwxr-xr-x 5 root root 124K Mar 15 11:02 ..

-rw-r--r-- 1 root root 9.6M Mar 15 11:01 endo.1.fa

-rw-r--r-- 1 root root 9.6M Mar 15 11:01 endo.2.fa

-rw-r--r-- 1 root root  80K Mar 15 11:01 segsites

ref.faと3つのサブディレクトリ(bact、cont、endo)ができる。bactはバクテリアコンタミ、contは同じ種からのコンタミ、endoは内性(endogenous)。endogenousのendo.1.faとendo.2.faのheterozygousポジションは以下のファイルにまとめられている。

> head data/endo/segsites 

# head data/endo/segsites 

>ref_1 312 C T

>ref_1 452 G T

>ref_1 1135 A G

>ref_1 1478 T A

>ref_1 1522 G T

>ref_1 2016 G C

>ref_1 2227 C T

>ref_1 4486 A C

>ref_1 5080 C G

>ref_1 5820 G A

 

2、リファレンスが用意できたらaDNAをシミュレートする。コマンドの最後には、1で作ったディレクトリを指定する。

cd ..

perl gargammel.pl -c 3 --comp 0,0.08,0.92 -f src/sizefreq.size.gz -matfile src/matrices/single- -o data/simulation data/

> ls -alh data/

# ls -alh data/

total 200M

drwxr-xr-x 1 root root 4.0K Mar 15 17:15 .

drwxr-xr-x 1 root root 4.0K Mar 15 10:58 ..

drwxr-xr-x 2 root root 4.0K Mar 15 11:01 bact

drwxr-xr-x 1 root root 4.0K Mar 15 17:13 cont

drwxr-xr-x 1 root root 4.0K Mar 15 17:13 endo

-rw-r--r-- 1 root root 9.6M Mar 15 11:01 ref.fa

-rw-r--r-- 1 root root    0 Mar 15 17:13 simulation.b.fa.gz

-rw-r--r-- 1 root root 1.3M Mar 15 17:13 simulation.c.fa.gz

-rw-r--r-- 1 root root  28M Mar 15 17:13 simulation.e.fa.gz

-rw-r--r-- 1 root root  36M Mar 15 17:14 simulation_a.fa.gz

-rw-r--r-- 1 root root  30M Mar 15 17:13 simulation_d.fa.gz

-rw-r--r-- 1 root root  48M Mar 15 17:14 simulation_s1.fq.gz

-rw-r--r-- 1 root root  49M Mar 15 17:14 simulation_s2.fq.gz

roo

simulation_s1.fq.gzとsimulation_s2.fq.gzが最終出力。

 

実行例

Githubには他にもいくつかの例がある。

 

1、非常に低いcoverageの40-bpリードを発生させる。

gargammel.pl -c 0.5 --comp 0,0,1 -l 40 -o data/simulation data/
  • -c     Endogenous coverage
  • --comp [B,C,E] Composition of the final set in fraction the 3 numbers represent the bacterial, contaminant and endogenous
    ex: --comp 0.6,0.02,0.38 will result in 60% bacterial contamination while the rest will be from the same species 5% will be contamination and 95% will be endogenous Default: --comp 0,0,1
  • -l     Generate fragments of fixed length (default: 35)
  • -o    Output prefix (default: [input dir]/simadna)

 

2、20xの高いcoverageだが現代のコンタミネーションが激しい(40%)45-bpのリードを発生させる。

gargammel.pl -c 20 --comp 0,0.4,0.6 -l 45 -o data/simulation data/

 

3、Briggs et al. 2007 model(pubmed)を使い、脱アミノ化した二本鎖DNAを1M発生させる。

gargammel.pl -n 1000000 --comp 0,0,1 -l 40 -damage 0.03,0.4,0.01,0.3 -o data/simulation data/
  • -n    Generate [number] fragments (default: 1000)
  • -damage [,,,s]   For the Briggs et al. 2007 model. The parameters must be comma-separated e.g.: -damage

 

他の例はGithubを確認して下さい。丁寧に説明されています。 

引用

gargammel: simulations of ancient DNA datasets

Gabriel Renaud, Kristian Hanghøj, Eske Willerslev, Ludovic Orlando
Bioinformatics, Volume 33, Issue 4, 15 February 2017, Pages 577–579