macでインフォマティクス

macでインフォマティクス

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

RepeatScoutでトランスポゾンなどのリピートをde novoで探す

 RepeatScoutはゲノム中のトランスポゾンなどのリピートを探すツール。リピートを見つけると、そのシードを保存性がなくなるまで伸長する戦略をとることで、見つかりにくい長くてやや配列に違いがあるリピートまで探索することが可能とされる(タンデムリピートやlow-complexityリピートは本手法のターゲットではない)。

 オーサーらの用意したデータでは、競合のRECONと比較して10倍以上短い時間で、2倍以上のリピートが検出されている。 RepeatMaskerの出力に対応している。

 

 

インストール

依存

  • Tandem Repeats Finder

https://tandem.bu.edu/trf/trf409.macosx.download.html

  • RepeatMasker

http://www.repeatmasker.org/RMDownload.html

Tandem Repeat Finderはバイナリをダウンロードして、下のようにリネームする。

 

Github

https://github.com/mmcco/RepeatScout

brewで導入できるが、サブコマンドが入らないので自分でビルドする。

git clone https://github.com/mmcco/RepeatScout.git
cd RepeatScout
make

フォルダ全体にパスを通しておく。Tandem Repeats Finderもここにコピーして、trfという名前にリネームする。

 

リピートライブラリ

http://bix.ucsd.edu/repeatscout/

 

ラン

ランは複数段階で行う。

1、データベースのビルド。全ての1-merの配列をpick upしてテーブルにする。

build_lmer_table -l 14 -sequence input.fasta -freq output.freq

2、そのテーブルファイルからFASTAを作る。

RepeatScout -sequence input.fasta -output output_repeats.fasta -freq output.freq -l 14

3、単純リピートなどを除外する。またデフォルトでは繰り返し数が10以下のリピートも排除する。

cat output_repeats.fasta | filter-stage-1.prl > repeats_filtered_stg1.fasta

4、RepeatMaskerでフィルタリングされた領域を分析する。

RepeatMasker -pa 20 -s -lib repeats_filtered_stg1.fasta input.fasta &

5、step4と並行して、規定回数登場しなかったリピートを排除する作業を行う。

cat repeats_filtered_stg1.fasta | filter-stage-2.prl --cat=Final_assembly.fasta.out --thresh=3 > repeats_filtered_stg2.fasta

6、RepeatMaskerで検出された部位を元に、step5の結果から最終的なリピート情報を出力する。

RepeatMasker -pa 20 -s -lib repeats_filtered_stg2.fasta input.fasta

 

 

 

 

 

 

引用

De novo identification of repeat families in large genomes.

Price AL1, Jones NC, Pevzner PA.

Bioinformatics. 2005 Jun;21 Suppl 1:i351-8.

 

SEQanswers

http://seqanswers.com/forums/showthread.php?t=5448

 

マイクロサテライトをraw readsから直接探すpalfinder

 palfinderはマイクロサテライトやsimple sequence repeats (SSRs)を探すツール。454やilluminaのNGSデータから直接マイクロサテライトを検出し、さらに内部でprimer3を動かし、その増幅プライマーを設計する機能を備える。

 

 

インストール

依存

  • primer3

primer3はbrewでインストールできる。

brew install primer3

 本体のダウンロード

https://sourceforge.net/projects/palfinder/?source=typ_redirect

 

ラン

ランはconfig.txtのパラメータを元に行われる。また、単純リピートをprimer3設計サイトから排除するために、単純リピートをリスト化したファイルも必要となる (simple.ref)。

configファイルは以下のようなものを記載する。

neo$ cat config.txt

 

# config.txt

 

findPrimers 1 #プライマーも設計したいなら1、マイクロサテライトのstatisticsだけなら0

platform Illumina #454かilljuminaを記載

inputFormat fastq

pairedEnd  1 #ペアリードなら1、そうでないなら0。454はシングル、illuminaはペアのみ対応

input454reads  test/data/454_All_python.fna #454データのパス。formatはFASTA(illuminaならこの項目は不要)

inputReadFile  test/data/IlluminaPE_SG_10k_1.fq #illuminaデータのパス。フォワードリード。scarf format。(454ならこの項目は不要)

pairedReadFile  test/data/IlluminaPE_SG_10k_2.fq #リバースリード

MicrosatSumOut  test/output/test_microsat_summary.txt #出力パス

PALsummaryOut  test/output/test_PAL_summary.txt #出力パス

2merMinReps 6 #2mer単位のリピートの最低繰り返し数。以下だと検出しない。

3merMinReps 0 #3mer

4merMinReps 0

5merMinReps 0

6merMinReps 0

primer3input  test/output/pr3in.txt #primer3一時ファイル1

primer3output  test/output/pr3out.txt #primer3一時ファイル2

keepPrimer3files  0 #primer3一時ファイル3

primer3executable  /usr/local/bin/primer3_core #プライマー3コアのパス

prNamePrefix  test_ #primer名。これだとtest_という名前になる。

 

 

#-------------------- primer3 Parameters ---------------------#

PRIMER_TASK  pick_pcr_primers

PRIMER_OPT_SIZE  20

PRIMER_MIN_SIZE  18

PRIMER_MAX_SIZE  30

PRIMER_MAX_NS_ACCEPTED  0

pr3ProductSizeRangeMinVal  60

pr3ProductSizeRangeMaxVal  500

PRIMER_OPT_SIZE  20

PRIMER_MIN_GC  30.0

PRIMER_MAX_GC  80.0

PRIMER_GC_CLAMP  2

PRIMER_MAX_END_GC  5

PRIMER_MIN_TM  58.0

PRIMER_MAX_TM  65.0

PRIMER_OPT_TM  62.0

PRIMER_PAIR_MAX_DIFF_TM  2.0

PRIMER_TM_FORMULA  0

PRIMER_MAX_SELF_ANY  8.00

PRIMER_PAIR_MAX_COMPL_ANY  8.00

PRIMER_MAX_SELF_END  3.00

PRIMER_PAIR_MAX_COMPL_END  3.00

PRIMER_MAX_POLY_X  4

PRIMER_LOWERCASE_MASKING  0 

PRIMER_NUM_RETURN  1

PRIMER_MISPRIMING_LIBRARY  simple.ref

PRIMER_MAX_LIBRARY_MISPRIMING  10.00

PRIMER_LIB_AMBIGUITY_CODES_CONSENSUS   0

 

ランはconfigファイルを指定して行う。 

perl pal_finder_v0.02.04.pl config.txt

 

出力

454リードからの探索結果。

f:id:kazumaxneo:20170925003448j:plain

マイクロサテライト予測結果サマリー。

f:id:kazumaxneo:20170925003607j:plain

 

  

引用

Rapid Microsatellite Identification from Illumina Paired-End Genomic Sequencing in Two Birds and a Snake

Todd A. Castoe,1 Alexander W. Poole,1 A. P. Jason de Koning,1 Kenneth L. Jones,1 Diana F. Tomback,2 Sara J. Oyler-McCance,3 Jennifer A. Fike,3 Stacey L. Lance,4 Jeffrey W. Streicher,5 Eric N. Smith,5 and David D. Pollock1,*

PLoS One. 2012; 7(2): e30953.

 

mrepsでタンデムリピートを探す

mrepsはダイレクトリピートを探すツール。短い単位の繰り返し配列がタンデムに続く領域を検出することができる。

 

ミニチュートリアル

http://mreps.univ-mlv.fr/tutorial.html

webサーバー版

http://bioinfo.lifl.fr/mreps/mreps.php

 

インストール

Github

GitHub - gregorykucherov/mreps: mreps: software for tandem repeat identification in DNA

git clone https://github.com/gregorykucherov/mreps.git
cd mreps
make
mreps -h

 

ラン

デフォルト条件のラン。

mpreps -fasta input.fa > output

 出力

user$ head -20 out 

 

 *****************************************************************************

 *                              mreps 2.6                                    *

 *                                                                           *

 *                Finding tandem repeats in DNA sequences                    *

 *                                                                           *

 *                      http://mreps.univ-mlv.fr/                            *

 *****************************************************************************

 

Processing sequence 'chr'

 

* Processing window [1 : 2695903] *

 

   from   ->       to  :  size <per.> [exp.] err-rate sequence

 ---------------------------------------------------------------------------------------------

    2111  ->      2122 :    12  <3>  [4.00]  0.000 GAC GAC GAC GAC 

   18777  ->     18795 :    19  <9>  [2.11]  0.000 GACAACCGC GACAACCGC G 

   20562  ->     20579 :    18  <9>  [2.00]  0.000 TGGCAGCAA TGGCAGCAA 

   21928  ->     21944 :    17  <6>  [2.83]  0.000 CGATCG CGATCG CGATC 

   30984  ->     30998 :    15  <6>  [2.50]  0.000 GCGATC GCGATC GCG

バクテリアゲノムなら数秒以内に解析できる。GACx4のリピートなどが検出されている。末端のリピートは不完全でも検出されるため、expの列は必ずしも整数倍にはならない。

 

不完全なリピートを探すこともできる。

1塩基の間違いを許容してダイレクトリピートを探す。

mpreps -res 1 -fasta input.fa > output

出力

usr$ head -20 out 

 

 *****************************************************************************

 *                              mreps 2.6                                    *

 *                                                                           *

 *                Finding tandem repeats in DNA sequences                    *

 *                                                                           *

 *                      http://mreps.univ-mlv.fr/                            *

 *****************************************************************************

 

Processing sequence 'chr'

 

* Processing window [1 : 2695903] *

 

   from   ->       to  :  size <per.> [exp.] err-rate sequence

 ---------------------------------------------------------------------------------------------

    2111  ->      2122 :    12  <3>  [4.00]  0.000 GAC GAC GAC GAC 

    3048  ->      3065 :    18  <6>  [3.00]  0.167 GAGCTG GATCTG GACCTG 

    3904  ->      3925 :    22  <9>  [2.44]  0.077 TGGTCGGCA TGGTCGGCT TGGT 

    4611  ->      4626 :    16  <7>  [2.29]  0.111 CCGCTGA CCGCCGA CC 

    8829  ->      8847 :    19  <6>  [3.17]  0.154 CGTTTG CGATTG CGATCG C

err-rateが0でないリピートも検出されている。

 

 

領域1000-12000まで、最低3回以上続くダイレクトリピートを探し、lessに渡す。

mreps -res 3 -exp 3.0 -from 10000 -to 12000 input.fa |less
  •  -s specifies the sequence in command line
  • -fasta allows DNA sequences in FASTA format
  • -res "resolution" (error level)
  • -from starting position n
  • -to end position n
  • -minsize repeats whose size is at least n
  • -maxsize repeats whose size is at most n
  • -minperiod repeats whose period is at least n
  • -maxperiod repeats whose period is at most n
  • -exp repeats whose exponent is at least x
  • -allowsmall output small repeats that can occur randomly
  • -win process by sliding windows of size 2*n overlaping

 

 

 

引用

mreps: efficient and flexible detection of tandem repeats in DNA

R. Kolpakov, G. Bana, and G. Kucherov

Nucleic Acid Research, 31 (13), July 1 2003, pp 3672-3678.

 

Finding approximate repetitions under Hamming distance

R. Kolpakov, G. Kucherov

Theoretical Computer Science, 2003, vol 303 (1), pp 135-156.

 

ShortStackでsmall RNAをアノテートする

ShortStackはsmall RNA seqのデータをリファレンスゲノムにアライメントし、small RNAのlociをアノテートするツール。改良が続けられており、2報目の論文では、高速化の他、複数のシーケンスデータの入力、bowtieによるアライメントなどに対応した。

 

 

テストデータ

https://psu.app.box.com/v/axtelldata/

 

インストール

依存

  • samtools (version 1.x or higher)
  • bowtie (if aligning)
  • bowtie-build (if aligning and .ebwt indices not found)
  • gzip (if aligning)
  • RNAfold (unless running with --nohp option to disable MIRNA search)

gzipmacに標準でインストールされている。他はbrewで導入できる。

本体

Github

https://github.com/MikeAxtell/ShortStack

バイナリーをダウンロードしてパスを通す。

git clone https://github.com/MikeAxtell/ShortStack.git
ShortStack -h #動作確認

 

ラン

ランにはsmall RNA seqのシーケンスデータ(fastq、fasta、color-space)が必要である。ペアリードには対応していないが、複数データがある場合、コンマで区切って入力可能。また、シーケンスデータはgz圧縮されていても使える。

 ShortStack --readfile input.fastq --genomefile reference.fa
  •  --genomefile path to reference genome in .fasta or .fa format. Required.
  • --outdir name of output directory to be created for results. Defaults to 'ShortStack_[time]', 
  • --readfile path to readfile(s) to be aligned. valid formats: .fasta, .fa, .fasta.gz, .fa.gz, .fastq, .fq, .fastq.gz, .fq.gz, .csfasta, .csfasta.gz. Multiple files, can be specified as separate arguments to --readfile ... e.g. --readfile file1.fastq file2.fastq file3.fastq Mutually exclusive with --bamfile or --cramfile.
  • --bamfile path to input .bam alignment file of small RNAs.
  • --cramfile path to input .cram alignment file of small RNAs.

アライメントが終わっている場合、fastqの代わりにbam (cram formatも可能)を指定することもできる。

出力

usr$ head -5 ShortStack_1506082238/Results.txt 

#Locus Name Length Reads RPM UniqueReads FracTop Strand MajorRNA MajorRNAReads Complexity DicerCall MIRNA PhaseScorShort Long 20 21 22 23 24

chr:4-52647 Cluster_1 52644 3099 11838.456 3076 0.504 . AACAGACCCUGAAAAUCCCAACUUCUCCAUUCCAUCCGGAGAGCAAAGAAGUAAGGGGGUUGAAUUCGAUAUCGCGGGGGAAAUCCUACCGGGCUGGAAUAUUAUUGCUUCCUAUGCUUAUACCGAUGCCAGGGUCACCAAGGAUGACAAUCUGGAGCCUGGUAAUUUGCUUGAGGGGGUUCCCUUUAACUCGGCCAGUUUGUGGUCAACUUACGAAAUUCAAGCCGGUGAUUUACAGGGUUUGGGAUUUGGCCUGGGAUUGUUUUAUGUGGGGGAACGCCAAGGUGAUUUAAAUAAUUCU 2 0.987 N N2 -1 0 3099 0 00

chr:52803-123897 Cluster_2 71095 4159 15887.751 4139 0.493 . GAUCGCCCCUUGGCCAGGGGAAUUCUCCUCCAGUGCUUGCAAGGGAGGGGCAAUAUAGGAAAAUACAAUCAACUCGAUCGCCGUCGAGCCGAAGUCGAGUAAAAACCGCUAUCAGGAGCCUCUAUGUACAUCGUUCAAAUUGCCUCAGAAUGCGCCCCCGUCAUUAAGGCUGGGGGAUUGGGGGAUGUUAUUUACGGCCUAAGCCGUGAAUUGGAACUGCGGGGCCAUUGCGUCGAGCUAAUCCUACCCAUGUACGAUUGCAUGCGCUAUGACCACAUCUGGGGUUUACACGAUGCUUACC 3 0.990 N N2 -1 0 4159 00

chr:124169-527481 Cluster_3 403313 24207 92472.896 24138 0.495 . AGGAAGCCAUUGAUCUGAUUAUUAAUGGCAUGCCGGUGCGGAGUAACUUAGAGUCAAAACUGUUCGGCAGCCAUACCCUUUCCUUGGCGAAAUCUACCAAAGUGCCGGUGAUGAUUUUACGCCCCCAAUUGGUCAGCACUUACACCGUUGAAGAAAUGGCUUUGCGGUGCCAACAUCUCUGGCGCAAUUUACUAGUGCCCUACGAUGCUAGUUCUGCGGGUAAUUAUUUAAUAGAAAGAUUAAAAAGUGCCUUGGAAAAGGCUCCCCCCGGUAAGGUUGAGUCCUGUUACUUCCUCUCCAU 3 0.989 N N2 -1 0 24207 00

chr:527956-1200296 Cluster_4 672341 40473 154610.466 40268 0.497 . AUACCCAUCCCCUAUGUUCAAUGGUCGGGGAACUGGUCCAAAUUGGCGAUCGUCUCUCAUCCGAUUUCCCUAAGUACACCAACUGCUCAAUAUCUGAAGCUUUAAAAUCAGGUUUUUAUAAAGUGGCGAUCGUUUGCUCAAUCUACUUUUGUAAAUCCUGUUCAUAAAGAGUUUUAGCCAUUAUAACUACCUUCUAUAUCAGGGCAAUUCAAUCCCAGUUGCUUUAGCAUAGCGUUUAGUGUUCCAGUUUUUAGUUCAUCACUAGGAUUUCGCACAAUGGUCAGGCGAUCGCUUCUACACU 3 0.989 N N2 -1 0 40473 00

 

 

 

 

Methodsの論文になっています。詳細は論文を確認してみて下さい。

 

引用

Identification and annotation of small RNA genes using ShortStack.

Shahid S1, Axtell MJ2.

Methods. 2014 May 1;67(1):20-7. doi: 10.1016/j.ymeth.2013.10.004. Epub 2013 Oct 16.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3989477/

 

メタゲノム向けの全ORF検出ツール OrfM

 

OrfMはcontigやアセンブルされていないリードからstopコドンの有無に関わらずorfを探索するツール。データサイズが莫大になるメタゲノム向けに設計された。非常に高速に動作し、translateやembossパッケージのgetorf、prodigalなどより数倍速く動作するとされる。例えばバクテリアゲノムのサイズなら1秒程度で全てのorfを探索する。エラーを考慮するような機能はないため、イルミナの高精度なシーケンスデータの使用が推奨されている。

 

 

インストール

Github

https://github.com/wwood/OrfM/releases

orfm-0.7.1.tar.gzをダウンロードして解凍する。

cd orfm-x.x.x 
./configure
make
sudo gem install rspec bio-commandeer
make check
sudo make install
orfm -h #ヘルプが出るか確認

 

ラン

デフォルトでは96-bp(-m 96)以上のorfを全て探す。stopやstartがなくても96bp以上は全て検出される(mは3の倍数の必要がある)。

orfm input.fasta >orfs.fa

 入力はfastaまたはfastqに対応している。

出力はアミノ酸配列となる。

>chr_2_5_1

TIHRGHIPPQIRLIHHVIVDQTSGVDHFGNFRQPAMAR

>chr_51_6_2

PNTKKIGQRGNPPLTLPPDCQINYSPWAHSAADPPHPPRHRGPN

>chr_1_4_3

VSKSVSVSKRLMLCKISPSAKSRITSDNNFTTSKLSKLAKYQEDWAKRKSPANTATRLSNKLFTVGTFRRRSASSTTSSWTKLAVWIISVISASRRWRA

>chr_293_5_4

PQMFTIFVKMITHIPGRIAQNPRCLGTVGNERGGFKIAFFQIGF

>chr_336_6_5

NLKYRPQSGTLTPNVYHIRKDDHTHPGAHRSKPPLPWHCWQ

ヘッダーのchrは使用したfastaのヘッダー名、次の数値はorfスタート位置、次はフレームの番号、最後は通し番号。

全てのフレームをチェックするのでポジションの重複も相当数出る。

 

  • -m minimum number of nucleotides (not amino acids) to call an ORF on [default: 96]
  • -t output nucleotide sequences of transcripts to this path [default: none]
  • -l ignore the sequence of the read beyond this, useful when comparing reads from with different read lengths [default: none]
  • -c ID codon table for translation (see http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes for details) [default: 1]
  • -p print the actual stop codons at sequence ends if encoded [default: do not]
  • -s only print those ORFs in the same frame as a stop codon [default: off]

 

 

引用

OrfM: a fast open reading frame predictor for metagenomic data

Ben J. Woodcroft Joel A. Boyd Gene W. Tyson

Bioinformatics, Volume 32, Issue 17, 1 September 2016, Pages 2702–2703,

NCBIからバクテリアゲノムをダウンロードする

 

コンプリートなゲノムのダウンロード。

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/assembly_summary.txt 
awk -F '\t' '{if($12=="Complete Genome") print $20}' assembly_summary.txt > assembly_summary_complete_genomes.txt
mkdir bacteria_genome
for next in $(cat assembly_summary_complete_genomes.txt); do wget -P bacteria_genome "$next"/*genomic.fna.gz; done
gunzip bacteria_genome/*.gz

 

 blastデータベースの作成

cat *.fna > all_complete_bacteria.fna
makeblastdb -in all_complete_bacteria.fna -parse_seqids -dbtype nucl -title bacteria -out bacteria
makeblastdb -in all_complete_bacteria.fna -parse_seqids -dbtype nucl -title bacteria -out bacteria
  •  -parse_seqids Option to parse seqid for FASTA input if set, for all other input types seqids are parsed automatically
  • -title Title for BLAST database Default = input file name provided to -in argument

 

 

Ensemblから全ゲノムをダウンロードする。

rsync -av rsync://ftp.ensembl.org/ensembl/pub/current_fasta/*/dna/*.dna.toplevel.fa.gz ./ 

 

 

 

引用

Biostar

https://www.biostars.org/p/61081/

Download All The Bacterial Genomes From Ncbi

 

https://nsaunders.wordpress.com/2013/05/28/how-to-bulk-retrieval-of-archaeal-genome-sequences-from-the-ncbi-ftp-site/

bamに塩基置換やindel変異を起こすbamsurgeon

 

bamsurgeonはガンの原因となる体細胞突然変異をシミュレートするために構築されたbamに対する変異導入ツール。ユーザーが用意したリストを元にして、bamに不完全な変異や構造変化を引き起こす大きな変異を導入することができる。2015年にnature methodsに発表された。

 

マニュアル

https://hpc.nih.gov/apps/bamsurgeon/Manual.pdf

インストール

Github

https://github.com/adamewing/bamsurgeon

git clone https://github.com/adamewing/bamsurgeon.git
python setup.py build
sudo python setup.py install

 

 

ラン

塩基置換の導入

SNV変異をbamに導入する (test/)。

addsnv.py -v random_snvs.txt -f testregion_bwamem_realign.bam -r Homo_sapiens_chr22_assembly19.fasta -o change.bam

 #リストについてはもう少し下の方に記載。

 

 

手術したbam出力を確認する。

igv -g Homo_sapiens_chr22_assembly19.fasta testregion_bwamem_realign.bam,change.bam

22 33714964 33714965 0.5  50%変異のはず(下が改変後)。

f:id:kazumaxneo:20170917201850j:plain

  • -v Target regions to try and add a SNV, as BED
  •  --mindepth minimum read depth to make mutation (default = 10)
  • --maxdepth maximum read depth to make mutation (default = 2000)
  • -m allelic fraction at which to make SNVs (default = 0.5)
  • -f sam/bam file from which to obtain reads
  • -r reference genome, fasta indexed with bwa index -a stdsw _and_ samtools faidx
  • -o .bam file name for output
  • -s maximum allowable linked SNP MAF (for avoiding haplotypes) (default = 1)
  • --ignoresnp  make mutations even if there are non-reference alleles sharing the relevant reads
  • --ignoreref make mutations even if the mutation is back to the reference allele
  • --force force mutation to happen regardless of nearby SNP or low coverage
  • --single input BAM is simgle-ended (default is paired-end)
  • --maxopen  maximum number of open files during merge (default 1000)
  • --requirepaired skip mutations if unpaired reads are present
  • -d allow difference in input and output coverage (default=0.9)
  • -n maximum number of mutations to try (default: entire input)

#リストについて

SNVのリストは以下のようなフォーマットで準備する。1-3列目はBEDフォーマットと同じである。4列目はなくても動くが、0-1の数値を入れることで変異の割合を指定できる。4列目がなければ50%となる。

user$ head random_snvs.txt 

22 34166720 34166721 1.0

22 33908770 33908771 0.15

22 33714964 33714965 0.5

22 33769483 33769484 0.75

22 33958087 33958088 0.01

22 34264774 34264775 0.25

22 33702084 33702085 0.2

22 34141184 34141185

22 33926586 33926587

BEDOPSを使えば、VCFから作ることができる。

フォーマット変換 VCF, GTF, GFF => BED - macでインフォマティクス

vcf2bed --snvs < input.vcf > SNV.bed
cut -f 1,2,3,4,6 SNV.bed |sed -e 's/\./0\.5/g' > SNV.txt

出力

user$ head -5 SNV.txt 

Chromosome 198597 198598 0.5 C

Chromosome 387005 387006 0.5 C

Chromosome 528665 528666 0.5 G

Chromosome 609078 609079 0.5 T

Chromosome 609083 609084 0.5 A

あとは0.5を必要な変異率に変えればリストに使うことができる。 

 

 

 

 

 small indelの導入 

small indel変異をbamに導入する(test/)。

addindel.py -v test_indels.txt -f testregion_bwamem_realign.bam -r Homo_sapiens_chr22_assembly19.fasta -o human_indel.bam

 SNVのリストは以下のようなフォーマットで準備する。indel変異では4列目以降も必ず必要である。

user$ cat test_data/test_indels.txt 

22 33694303 33694304 0.5 INS ATGGC

22 33815273 33815274 0.25 INS CG

22 33920490 33920491 0.33 INS ATTTGCTTAGCTGAGGGCTTAGGCTTAGCATGCGTA

22 33944542 33944543 0.50 DEL

22 34006643 34006653 0.50 DEL

22 33958087 33958090 0.40 DEL

5列目はINSかDELのみ許されている。挿入は6列目に挿入される配列を記載する。

 

やはりVCFからリストを作ることができる。

#insertion
vcf2bed --insertions < input.vcf > ins.bed #insertionのみbedに変換
cut -f 1,2,3,4 ins.bed |sed -e 's/\./0\.5/g' > temp
cut -f 4,6 ins.bed |sed -e 's/\./INS/g' |paste temp - > insertion.txt

#deletion
vcf2bed --deletions < input.vcf > del.bed
cut -f 1,2,3,4 del.bed |sed -e 's/\./0\.5/g' > temp
cut -f 4,6 del.bed |sed -e 's/\./DEL/g' |paste temp - > deletion.txt

#マージしてソートする。
cat insertion.txt >> deletion.txt
sort -k 1,1 -k 2n deletion.txt > indel.txt

 

 

del

22 33944542 33944543 0.50 DEL

f:id:kazumaxneo:20170917233108j:plain

 

 In

f:id:kazumaxneo:20170917222029j:plain

 igvで見るには、ソートしてindexをつけておく必要がある(e.g.,  samtools sort i.bam > o.bam && samtools index o.bam) 

 

 

 

SVの導入 

構造変化(SV)を引き起こす大きな変異も導入できる(insertions, deletions, duplications, inversions, and compound rearrangements)。

addsv.py -p 18 -v test_sv2.txt -f testregion_bwamem_realign.bam -r Homo_sapiens_chr22_assembly19.fasta -o human_sv.bam 

 

SVのリストは以下のようなフォーマットで準備する。

user$ cat test_data/test_sv.txt 

22 34171055 34184700 INS ACTAACCTGCACAATGTGCACATGTACCCTAAAACTTAGAGTATAATAAAAAAAA 13 AA^TTTT

22 33607508 33607508 TRN 22 34314981 34314981

22 33871043 33884754 DEL 0.9

22 34071043 34084754 INV

22 34271043 34284754 DUP 3

細部は異なるが、small indelと同じような方法でVCFからリストは準備できる。ただしSV検出ツールによっては、VCFに変異後の配列を記載しないものがある。大きな挿入でも塩基を記載するようなVCF(例えばpindelならO.K)を変換元にしてください。

 

 

large deletion。

22 33871043 33884754 DEL 0.9

f:id:kazumaxneo:20170918001112j:plain

指定領域をアセンブルして、できた最長のcontigの9割(0.9)の領域が欠損となる。 

 

duplication。

22 34271043 34284754 DUP 3

f:id:kazumaxneo:20170918001838j:plain

3回複製される。アセンブルしているので、duplicationが起きる部位は、指定領域より短くなる(上の図)。

 

逆位

22 34071043 34084754 INV

f:id:kazumaxneo:20170918002007j:plain

上流側境界。

 

ACTAACCTGCACAATGTGCACATGTACCCTAAAACTTAGAGTATAATAAAAAAAAを挿入。ただし13-bp末端を複製して挿入する。また、AA^TTTTが領域内になれば、そこをターゲットとする(AAとTTTTの間に導入)

22 34171055 34184700 INS ACTAACCTGCACAATGTGCACATGTACCCTAAAACTTAGAGTATAATAAAAAAAA 13 AA^TTTT 

f:id:kazumaxneo:20170918004314j:plain

13-bp末端の複製があるので、13-bpズレて挿入部位が見えている。ターゲット配列を書かなければ、アセンブルされた最長のcontigの真ん中に導入される。

 

 

同時に複数の変異を導入、例えば欠損を起こしてそこに挿入を起こすような変異も導入できます。詳しくはマニュアルを参照してください。

https://hpc.nih.gov/apps/bamsurgeon/Manual.pdf

 

引用

Combining tumor genome simulation with crowdsourcing to benchmark somatic single-nucleotide-variant detection

Adam D Ewing, Kathleen E Houlahan, Yin Hu, Kyle Ellrott, Cristian Caloian, Takafumi N Yamaguchi, J Christopher Bare, Christine P'ng, Daryl Waggott, Veronica Y Sabelnykova, ICGC-TCGA DREAM Somatic Mutation Calling Challenge participants, Michael R Kellen, Thea C Norman, David Haussler, Stephen H Friend, Gustavo Stolovitzky, Adam A Margolin, Joshua M Stuart & Paul C Boutros

Nature Methods 12, 623–630 (2015) doi:10.1038/nmeth.3407