macでインフォマティクス

macでインフォマティクス

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

Complete Genomicsのシーケンスリードをマッピングする sirfast

 

 ハイスループットシークエンシング(HTS)技術は、[論文より ref.1]におけるペアエンド配列決定、および全ゲノムショットガンシーケンシング(WGS)[ref.2]の最初の使用以来、魅力的な速度で進化し続けている。 Roche / 454 [ref.3]、Illumina [ref.4]、ABI  SOLiD [ref.5]、Pacific Biosciences [ref.6]など多くのシーケンシング技術が開発されている[ref.7]。これらのプラットフォームの多くでは、リードのマッピング[ref.8]、バリアントのコール[ref.9,10]、ゲノムアセンブリ[ref.11]のための多くのツールが開発されている。

 多くの異なるアルゴリズムとそのオープンソース実装が利用可能であるため、研究者はHTS技術によって生成されたデータをより有効に利用することができる。しかし、Complete Genomics [ref.12](CG)プラットフォームによって生成されたデータをマップし分析するために考案された公に利用可能なアルゴリズムは存在しない。これは主に、Complete Genomics社がシーケンシングサービスのみを提供し、独自のソフトウェアを使用して独自の分析結果を顧客に提供しているためである。 CGデータを分析するために考案された唯一のアルゴリズムセットは、同社の科学者グループによって記述されている[ref.13]が、これらのアルゴリズムの実装は独自のものであり、未発表のままである。間違いなく、同社独自の分析パイプラインは、異なる特性を持つデータに対してきめ細かく調整されている。それでも、1000 Genomes Project [ref.14]のような多くのゲノムの解析と新しいアルゴリズム開発を組み込んだ研究プロジェクトは、解析結果をさらに改善するためのオープンソースのツールとアルゴリズムが利用できるというメリットがある。 CGデータを分析するためのオープンソースツールがあれば、他の研究者が他の方法で発見することができない新しい結果を開発することが可能になる。

 本論文では、CGプラットフォームで生成されたデータ用に設計された新しいリードマッパーsirFASTを紹介する。 Complete Genomicsは、他のすべてのHTS技術とは異なり、リードは複数のサブリード(読み取りセクション)で構成され、重複またはギャップを持つ可能性がある。著者らはsirFASTを設計して、リードごとにフレキシブルな "予想される"ギャップ/オーバーラップが存在する場合に、そのようなリードをリファレンスアセンブリに効率よくアライメントするにした。 Burrows-Wheeler変換[ref.15]に基づくメソッドはindelsでうまく拡張できないため、sirFASTをハッシュベースのシード・アンド・延長アルゴリズムとして実装した。イルミナとSOLiDをベースにしたツール[ref.16,17,18,19]と同様に、sirFASTはSAM [ref.20]とDIVET [ref.21]の両方のファイル形式をサポートしている。 

 

インストール 

Github

https://github.com/BilkentCompGen/sirfast

git clone https://github.com/BilkentCompGen/sirfast.git
cd sirfast/
make

./sirfast 

 $ ./sirfast 

sirFAST : Short interrupted Read Fast Alignment Search Tool. Enhanced with FastHASH.

 

Usage: sirfast [options]

 

General Options:

 -v|--version Current Version.

 -h Shows the help file.

 

 

Indexing Options:

 --index [file] Generate an index from the specified fasta file. 

 

 

Searching Options:

 --search [file] Search in the specified genome. Provide the path to the fasta file. 

Index file should be in the same directory.

 --pe Search will be done in Paired-End mode.

 --seq [file] Input sequences in fasta/fastq format [file]. If 

paired end reads are interleaved, use this option.

 --seq1 [file] Input sequences in fasta/fastq format [file] (First 

file). Use this option to indicate the first file of 

paired end reads. 

 --seq2 [file] Input sequences in fasta/fastq format [file] (Second 

file). Use this option to indicate the second file of 

paired end reads.  

 -o [file] Output of the mapped sequences. The default is "output".

 -u [file] Save unmapped sequences in fasta/fastq format.

 --best   Only the best mapping from all the possible mapping is returned.

 --seqcomp Indicates that the input sequences are compressed (gz).

 --outcomp Indicates that output file should be compressed (gz).

 -e [int] Maximum allowed hamming distance (default 4% of the read length).

 --min [int] Min distance allowed between a pair of end sequences.

 --max [int] Max distance allowed between a pair of end sequences.

 --maxoea [int] Max number of One End Anchored (OEA) returned for each read pair.

We recommend 100 or above for NovelSeq use. Default = 100.

 --maxdis [int] Max number of discordant map locations returned for each read pair.

We recommend 300 or above for VariationHunter use. Default = 300.

 --crop [int] Trim the reads to the given length.

 --sample [string] Sample name to be added to the SAM header (optional).

 --rg [string] Read group ID to be added to the SAM header (optional).

 --lib [string] Library name to be added to the SAM header (optional).

 

パスの通ったディレクトリに移動しておく。 

 

ラン

テストデータを使う。はじめにindexファイルを作成する。 

sirfast --index Sample_Data/sample_genome.fa 

 

シングルエンドのテストデータをsample_genomeにマッピングする。シーケンスデータはtsv形式で、以下のようになっている(format)。

f:id:kazumaxneo:20180413214026j:plain

マッピングを実行する。

sirfast --search Sample_Data/sample_genome.fa --seq Sample_Data/sample_reads_e0 -e 0 -o Sample_Data/out.sam
  • --seq    Input sequences in fasta/fastq format [file]. If --seq [file] Input sequences in fasta/fastq format [file]. If  paired end reads are interleaved, use this option.
  • -e    Maximum allowed hamming distance (default 4% of the read length).
  • -o    Output of the mapped sequences. The default is "output".

out.sam (mapped) とunmappedが出力される。デフォルトでは条件を満たした全ての部位にマッピングされる。ベストマッピングだけに限定したければ--bestをつける。

 

 

引用

Fast and Accurate Mapping of Complete Genomics Reads

Donghyuk Lee,†, Farhad Hormozdiari,†, Hongyi Xin,a Faraz Hach, Onur Mutlu, and Can Alkand

Methods. 2015 Jun 1; 0: 3–10.

 

What is Complete Genomics data and how does it compare with other genomics datasets?

What is Complete Genomics data and how does it compare with other genomics datasets? - Quora

germlineとsomaticの変異を検出する SNVSniffer

 

 次世代シークエンシング(NGS)に基づいて、単一ヌクレオチド変異(SNV)または挿入 - または欠失(indel)突然変異を呼び出すための様々なアプローチが開発されている。しかし、それらの大部分は、特定のタイプの突然変異に捧げられている。正常細胞における生殖系列SNV、癌/腫瘍細胞における体細胞SNV、またはindelsのみである。文献では、生殖系列および体細胞SNVs / indelsの両方に対する効率的かつ統合されたコーラーは、まだ広範に調査されていない。

 SNVSnifferは、NGSデータからの生殖系列と体細胞SNVs / indelsの両方を同定する効率的かつ統合されたコーラーである。このアルゴリズムでは、ベイジアン確率モデルを使用してSNVを特定し、indelsを呼び出すため複数のungappedアライメントアプローチを調べることを提案する(一部略)。

 多様なベンチマークを使用して、包括的なパフォーマンス評価を実施した。生殖系列変異コールのために、SNVSnifferは最先端のFaSD、GATKおよびSAMtoolsと比較して優れたスピードで競争力の高い正確性を示す。体細胞変異のコールのために、このアルゴリズムは、VarScan2、SomaticSniper、JointSNVMix2、MuTectよりも高速で、同等またはそれ以上の精度を達成する。

 

 

本ツールには、germlineとsomaticの変異を再現可能なシーケンスリードのシミュレータ機能も備わっています。ユーザー定義の条件で変異を発生させて、本ツールで解析し、評価する、という流れが本ツールのみで可能になっています。

 

公式サイト 

http://snvsniffer.sourceforge.net/homepage.htm#latest

 version2マニュアル

http://snvsniffer.sourceforge.net/version2.htm#latest

 

インストール

cent OSに導入した。

依存

  • SAMtools

 

本体は上記公式から実行ファイルをダウンロードして使う。SNVSniffer 1.0はSNVのみがターゲットだが、SNVSniffer 2.0 はindelも検出可能になっている。

$ ./SNVSniffer 

 

SNVSniffer (v2.0.4): identifying germline and somatic single nucleotide and indel mutations (released on March 07, 2016)

Command:

snp            Identify single-sample SNVs/indels

somatic        Identify germline/somatic SNVs/indels from tumor-normal pairs

gsim           Simulate sample reads with germline SNVs/indels

ssim           Simulate normal and tumor sample reads with somatic SNVs/indels

eval           Evaluate predicted germline/somatic SNVs/indels against gold standard, only for evaluation purpose

 

  

 

ラン

5つのコマンドのうち、performance結果を評価するevalコマンド以外の4つを紹介する。テストデータがアップされているので、それを使ってテストランする。リファレンスは"Home / data / v1.0 / synthetic_tumors_error_0.01"にある(リンク)。

 

1、snp - Identify single-sample SNVs/indels

$ SNVSniffer snp

 

Usage: SNVSniffer snp [options] sam_header infile

* sam_header: a SAM header file associated with the BAM file

* infile: a mpileup/pileup/BAM file (can have no header)

Input:

-g <string> (reference genome file, [Required by BAM format])

-f <int> (input file format, default = 2)

    0: mpileup format generated by SAMtools

    1: pipeup format generated by MAQ

    2: BAM file

 

Output:

-o <string> (output file name, default = STDOUT)

 

Base call and coverage:

-min_cov <int> (minimum coverage, default = 3)

-max_cov <int> (maximum coverage, default = 65536)

-min_bqual <int> (minimum Phred base quality score, default = 13)

-min_mapq <int> (minimum mapping quality score, default = 0)

-seq_err_rate <float> (sequencing error rate of the data, default = 0.01)

 

SNV calling:

-exe_mode <int> (execution model for SNV calling, default = 0)

    0: very fast and lesst accurate

    1: fast and accurate

    2: slow and most accurate

-prior <int> (model used for genotype prior probabilities, default = 1)

    0: equal probability

    1: prior probabilities with no consideration of Ti/Tv ratio

    2: prior probabilties by considering Ti/Tv ratio

-snp_rate <float> (SNP mutation rate for the species, default = 0.001)

-min_allele_freq <int> (minimum allele frequency, default = 0.2 [Important])

-stringency <int> (stringency level [0, 9] for low-confidence mutations, default = 6)

-homo_freq <float> (allelic frequency threshold for homozygous genotype, default = 0.75)

-relaxed <int> (using the relaxed confidence score calcuation (favours sensitivity), default = 1)

-min_locus_dist <int> (minimum distance interval between two neighboring SNP loci, default = 1 [0 disables it])

-use_strand_dist <bool> (use strand distribution information, default = 0)

-local_range <int> (set the local search range for strand distribution, default = 1000)

-call_snps <bool> (enable the calling of SNPs, default = 1)

-call_indel <bool> (enable the calling of insertions or deletions, default = 1)

-pvalue <float> (P value threshold for variant calling, default = 0.05)

ランにはSAMのヘッダー情報が必要。ヘッダーファイルを作成してからランする。

samtools view -H input.bam > header.sam 
SNVSniffer snp -f 2 header.sam -g reference.fa input.bam > output.vcf

出力。

$ head -n 15 out.vcf

##fileformat=VCFv4.1

##source=SNVSniffer-v2.0.4

##INFO=<ID=DP,Number=1,Type=Integer,Description="Coverage">

##INFO=<ID=HC,Number=1,Type=Float,Description="High confidence">

##INFO=<ID=LC,Number=1,Type=Float,Description="Low confidence">

##INFO=<ID=VT,Number=1,Type=String,Description="Variant type: SNP, INS or DEL">

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE

chr21 5010446 . T C . PASS DP=16;VT=SNP;HC=3.657 GT 0|1

chr21 5012778 . C G . PASS DP=26;VT=SNP;HC=4.350 GT 0|1

chr21 5013197 . C T . PASS DP=27;VT=SNP;HC=3.657 GT 0|1

chr21 5015438 . C A . PASS DP=45;VT=SNP;HC=4.350 GT 0|1

chr21 5016893 . C G . PASS DP=33;VT=SNP;HC=4.350 GT 0|1

chr21 5017812 . G A . PASS DP=31;VT=SNP;LC=3.657 GT 0|1

chr21 5018805 . G A . PASS DP=32;VT=SNP;HC=3.657 GT 0|1 

 

 

2、somatic - Identify germline/somatic SNVs/indels from tumor-normal pairs

$ SNVSniffer somatic 

 

Usage: SNVSniffer somatic [options] normal_sam_header tumor_sam_header normal tumor

1. normal_sam_header: SAM header file for normal

2. tumor_sam_header: SAM header file for tumor

3. normal: a mpileup/pileup/BAM file for normal

4. tumor: a mpileup/pileup/BAM file for tumor

Input:

-g <string> (reference genome file, required for BAM format input)

-f <int> (input file format, default = 2)

    0: mpileup format generated by SAMtools

    1: pipeup format generated by MAQ

    2: BAM format

 

Output:

-o <string> (output file name, default = STDOUT)

-o_somatic <int> (output somatic mutations, default = 1)

-o_loh <int> (output loss of heterzygosity mutations, default = 1)

-o_germline <int> (output germline mutations, default = 0)

-o_unknown <int> (output unknown type mutations, default = 1)

-o_indel <int> (output indels, default = 1)

 

Base call and coverage:

-min_cov_normal <int> (minimum normal coverage, default = 3)

-min_cov_tumor <int> (minimum tumor coverage, default = 3)

-min_bqual <int> (minimum Phred base quality score, default = 13)

-min_mapq <int> (minimum mapping quality score, default = 0)

-seq_err_rate <float> (sequencing error rate of the data, default = 0.01)

 

 

SNV calling:

-exe_mode <int> (execution model for SNV calling, default = 0)

    0: very fast

    1: fast and more accurate

    2: slow and most accurate

-prior <int> (model used for genotype prior probabilities, default = 1)

    0: equal probability

    1: prior probabilities with no consideration of Ti/Tv ratio

    2: prior probabilties by considering Ti/Tv ratio

-somatic_rate <float> (somatic mutation rate, default = 0.01)

-tumor_purity <float> (estimated purity (tumor data) of tumor sample, default = 0 [0 means AUTO])

-min_allele_freq <float> (minimum allele frequency for the normal, default = 0.2)

-stringency <int> (stringency level [0, 9] for low-confidence mutations, default = 6)

-homo_freq <float> (allelic frequency threshold for homozygous genotype, default = 0.75)

-p_value <float> (P value threshold for variant calling per sample, default = 0.05 [0 means global])

-use_strand_dist <bool> (use strand distribution information, default = 0)

-local_range <int> (set the local search range for strand distribution, default = 1000)

 

 

最初にnormalとtumorのbamからheader.samを作る。

samtools view -H chr21_e15_c30_normal.realigned.bam> normal_header.sam
samtools view -H chr21_e15_c30_tumor.realigned.bam > tumor_header.sam

somaticコマンドを走らせる。tumorの純度を"-tumor_purity <float>"で指定しない時はSNVSnifferによって推定された値が使われる。

SNVSniffer -f 2 -g chr21.fa somatic normal_header.sam tumor_header.sam chr21_e15_c30_normal.realigned.bam chr21_e15_c30_tumor.realigned.bam -o out.vcf

 

 

3、gsim - Simulate sample reads with germline SNVs/indels

$ SNVSniffer gsim

 

SNVSniffer gsim: a diaploid germline mutation simulator derived from WGSIM in SAMtools

Usage: SVNSniffer gsim [options] reference.fa reads.base

reference.fa: reference file stored in FASTA/FASTQ format

reads.base: the base file name for reads and VCF-formated mutations

Options:

Output:

-z (compress all of the output using ZLIB)

Random numbers:

-p <int> (random number seed value for generating the simulated genome, default=11)

-q <int> (random number seed value for generating simulated reads from the simulated genome, default = 47)

Reads:

-c <int> (coverage of the genome, default = 0 [>0 distables option -n])

-n <int> (number of read pairs, default = 1000000)

-1 <int> (length of the first read, default = 100)

-2 <int> (length of the second read, default = 100)

-a <float> (disgard if the fraction of ambiguous bases higher than #FLOAT, default = 0.05)

Germline:

-e <float> (sequencing base error rate, default = 0.02)

-i <int> (average insert size, i.e. outer distance, between the two ends, default = 500)

-d <int> (standard deviation of insert size, default = 50)

-r <float> (rate of mutations, including substitutions and indels, default = 0.001)

-R <float> (fraction of indels, default = 0.15)

-X <float> (probability an indel is extended, default = 0.3)

-H <float> (homozygous variant ratio, default = 0.3333)

-h (haplotype mode (all reads are sequenced from a single sequence)

カバレッジ30で、SNVとindelが含まれたgermlineのfastqをシミュレートする。

SNVSniffer gsim reference.fa output -e 0.01 -c 30

ランが終わるとoutput_1.fqとoutput_2.fqができる。下のキャプチャ画像は bamに変換してIGVで表示してみたもの。

f:id:kazumaxneo:20180412223627j:plain

 

 

4、ssim - Simulate normal and tumor sample reads with somatic SNVs/indels

$ SNVSniffer ssim

 

SNPSniffer ssim: a diaploid somatic mutation simulator

Usage: SNVSniffer ssim [options] reference.fa reads.base

reference.fa: reference file stored in FASTA/FASTQ format

reads.base: the base file name for reads and VCF-formated mutations

Options:

Output:

-z (compress all of the output using ZLIB)

Random numbers:

u <int> (random number seed value for the simulated genome, default=11)

v <int> (random number seed value for generating somatic mutations on the simulated gneome, default=149)

w <int> (random number seed value for generating reads from the simulated tumor genomes, default=97)

Reads:

-c <int> (coverage of the genome, default = 0 [>0 distables option -n])

-n <int> (number of read pairs, default = 1000000)

-1 <int> (length of the first read, default = 100)

-2 <int> (length of the second read, default = 100)

-a <float> (disgard if the fraction of ambiguous bases higher than #FLOAT, default = 0.05)

-i <int> (average insert size, i.e. outer distance, between the two ends, default = 500)

Somatic:

-E <float> (sequencing base error rate from simulated tumor genome, default = 0.02)

-s <float> (somatic mutation rate, default = 0.01)

-t <float> (ratio of homozygous variants (the rest are heterozygous), default = 0.3333)

-p <float> (probabiltiy of observing tumor reads at any somatic site [bionomial distribution], default = 0.9)

-q <float> (fraction of somatic indels, default=0.15)

-x <float> (probability an somatic indel is extended, default = 0.3)

-h <float> (homozygous rate for somatic variants, default = 0.3333)

Germline:

-e <float> (sequencing base error rate from simulated normal genome, default = 0.02)

-d <int> (standard deviation of insert size, default = 50)

-r <float> (rate of mutations, including substitutions and indels, default = 0.001)

-R <float> (fraction of indels, default = 0.15)

-X <float> (probability an indel is extended, default = 0.3)

-H <float> (homozygous variant rate, default = 0.3333)

 

 使い方はgsimと同様。

SNVSniffer ssim reference.fa reads_base -e 0.01 -c 30

 

 

 

 

参考

SNVSnifferの解説ページに、SNVSniffer以外のツールについても、normalとtumorを比較しtumor specificな変異を検出するコマンドが記載されている。覚書としてここに載せておきます。

準備 

正常組織。

#解説ページではbwa memでなくbwa alnでマッピングしている。
bwa aln -t 4 reference.fa normal_1.fa.gz >tmp1.sai
bwa aln -t 4 reference.fa normal_2.fq.gz >tmp2.sai

#saiからpaired.samを作り、bamに変換。
bwa sampe -r '@RG id:gatech SM:illumina' -a 700 reference.fa tmp1.sai tmp2.sai normal_1.fq.gz normal_2.fq.gz | samtools view -bhS -@ 10 - >normal.bam

#sortとindex
samtools sort -@ 10 normal.bam > normal_sorted.bam
samtools index normal_sorted.bam

#pileupは時間がかかる
samtools mpileup -s -C 50 -f reference.fa normal_sorted.bam >normal.pileup

腫瘍サンプルのマッピングとpileupも同様の手順で進める。

 

tumorの変異コール。

 Varscan

java -Xmx4g -jar VarScan.jar somatic normal.mpileup tumor.mpileup out.vcf -output-vcf 1 --strand-filter --tumor-purity 0.9 --min-coverage-normal 3 --min-coverage-tumor 3

JointSNVMix ( paper

java -Xmx4g -jar VarScan.jar somatic normal.mpileup tumor.mpileup out.vcf -output-vcf 1 --strand-filter --tumor-purity 0.9 --min-coverage-normal 3 --min-coverage-tumor 3

SomaticSniper (HP

bam-somaticsniper -F vcf -f reference.fa tumor.sorted.bam normal.sorted.bam out.vcf

Mutect

java -Xmx16g -jar muTect.jar --analysis_type MuTect --reference_sequence reference.fa --input_file:normal normal.sorted.bam --input_file:tumor tumor.sorted.bam --vcf out.vcf

SNVSniffer

SNVSniffer somatic normal_header.sam tumor_header.sam normal.sorted.pileup tumor.sorted.pileup -o out.vcf

 

引用

SNVSniffer: an integrated caller for germline and somatic single-nucleotide and indel mutations

Liu Y, Loewer M, Aluru S, Schmidt B

BMC Syst Biol. 2016 Aug 1;10 Suppl 2:47

 

リファレンスフリーで家族内変異や病変組織の変異を調べ、数十以下まで候補を絞り込む DIAMUND

追記 4/16

エラーが大量に出たので内容を修正しました。

 

   遺伝性疾患と癌の両方を含む、疾患の原因である突然変異を発見するためのゲノムシーケンシングの使用は、近年爆発的に増加している。全ゲノムシーケンスおよび全exomeシーケンスは、疾患表現型の原因となる可能性のある遺伝的変化を検出するために、何千人もの個体に採用されている。多くの成功は、ミラー症候群(Ng et al、2010)(わずか4人の脾臓の配列決定後に発見された)、バーター症候群[Choi et al、2009]の突然変異を同定するためのexomeシーケンシングを使用した2009年の報告に続き、フリーマン - シェルドン症候群[Ng et al、2009]に続いている。

   現在、高品質でほぼ完全なGRC37 / hg19ヒトリファレンスゲノム[The International Human Genome Sequencing Consortium、2001]は、これらの研究のすべてにおいて重要な役割を果たしている。分析のための標準プロトコールは、発端者および家族からのすべての配列をGRC37にアライメントさせ、次いで、一塩基多型(SNP)、小さな挿入および欠失(indels)、およびコピー数バリアントを検出する。しかし、ヒト個体群の自然変動のために、無作為の個体とGRC37とを比較すると、エクソームのみでも5万〜10万個の変異が生じる。ゲノム全体のシーケンスでは、数百万もの変異をもたらす。例えば、Roach et al(2010)は、ファミリーカルテットとヒトリファレンスゲノムを比較し450万個のSNPを見つけている。この膨大な候補の中で、しばしば単一の突然変異であるメンデリア病を引き起こす変異を見出すことは、乾草の中で針を探す (needle in a haystack) に類似したものになってしまう。最近の研究 (論文執筆時点) では、Strelka [Saunders et al、2012]およびFamSeq [Peng et al、2013]のような候補の膨大なセットを減らす統計モデルの開発に焦点が当てられてきたが、本質的な問題は残っている。リファレンスゲノムへのマッピングによって解析を始めると、非常に多数のリファレンスとの違いを検出することになる。

   Exome分析アルゴリズムでは、一連のフィルタを使用して50,000種類以上のバリアントのセットを管理可能なセットに絞り込み、個別に評価し、フォローアップ実験で検証することができる。しかしこれらのフィルタは、時には対象となる真のバリアントを除外する。例えば、フィルタリングの基準では、dbSNPのような大きなSNPデータベース、または1000Genomes Projectのような大きなプロジェクトで頻度> 1%で観察される変異すべてを除外する。これはフィルタとして非常に効果的であり、共通の母集団バリアントを表す何千もの無害なSNPを除去できる。しかし、このステップは、これらのデータベースのSNPのどれもが病気を引き起こすものではないことを暗黙のうちに仮定している。不完全な浸透、被験者の認識不能な疾患、およびその他の要因のために、健常者のみからSNPが収集されたとしても、この仮定は誤っている可能性がある。 解析パイプラインは、通常、多数の変異を除去するが、機能的に重要であるかもしれないスプライス部位および非コード変異も排除する。これらのフィルタや他のフィルタは意図せずに真の病因を排除する恐れがあり、さらにフィルタリング処理後でもexomeシーケンスで見つかる突然変異の数はまだ圧倒的に多いのが普通である。この問題は全ゲノムシークエンシングではさらに悪化し、無関係な変異の数は50倍になる。

   この論文では、exomeおよび全ゲノム解析に異なるアプローチをとるDiamund (direct alignment for mutation discovery) という新しい方法を紹介する。この方法では、突然変異候補は大幅に少なくなる。全てのサンプルをリファレンスゲノムにマッピングさせるのではなく、リードを互いに直接アライメントさせる。この方法は、(1)病変組織を同一個体の正常組織と比較する自己比較、および(2)被験体からのDNA配列間の差異がリファレンスとの比較よりずっと少ない家族内での解析で特に上手く機能する。 この方法では、germlineのリード (一般的にexome全体で1億以上の数)をGRC37リファレンスゲノムに合わせる必要はなく、複雑なゲノムアセンブリやこれらの大規模データセットからの統計的フィルタリングも使わない。以下で詳細に説明するように、著者らはどのサンプルでもユニークなシーケンスを素早く見つけることができる、より効率的なアルゴリズムを使用する。著者らはDiamundを2種類のexomes解析でテストした。まず、罹患した個体の病原組織に由来する初代培養線維芽細胞からのDNAを、同一個体の非病原性の初代培養線維芽細胞のDNAと比較したデータを検討した。腫瘍細胞または他の体細胞モザイクの遺伝的異常の分析のために、この直接的な比較は、最初に全てのリードをリファレンスゲノムと比較する分析よりも小さな変異数になる。第二に、family triosのデータを使った。このケースでは、子供のde novo突然変異が病気を引き起こしている疑いがあった。標準的なアルゴリズムでは、3つ全ての個体をリファレンスゲノムと比較するため、非常に大きな変異リストが出来上がるが、その多くは子供と親によって共有されるはずである。子供のDNAを両親と直接比較することにより、感度を失うことなく、またプロセスにノイズを加える家族特有の共通変異を検出することなく、すべてのde novo突然変異を迅速に同定することができた。真の変異を除いてしまう恐れのあるフィルターに頼る方法を取り除くことで、これらのテストでは、真のde novo変異リストは非常に少なかった。

   De novo突然変異は、メンデル症の高い割合を占める可能性がある。 Yang et alは250人の発端者(probands)およびその家族のexomeシークエンスについて調べ、常染色体優性患者33人およびX連鎖病の9人を同定した。これらのうち、常染色体優性の83%およびX連鎖変異の40%が新たに発生したものだった。ファミリー内のサンプル間、または影響を受けた組織と影響を受けていない組織との間の直接的な比較は、偽陽性を減らし、リファレンスゲノムから完全に欠けている領域の突然変異の検出も可能にする。いくつかのヒト集団は、多くのメガベースにまたがる大きなゲノム領域を共有していることが既に示されている[Li et al、2010]。これはヒトリファレンスゲノムから完全に欠けている。これらには、新規の断片的重複[Schuster et al、2010]ならびに全く新しい配列が含まれる。関心のある突然変異がこれらの領域の1つに落ちた場合、従来の方法はそれを見逃す。対照的に、著者らの直接比較アルゴリズムは、これらの領域も除かないため、それらの領域内での突然変異を見つけることができる。 重要な注意点として、Diamundは変種検出のよりgeneralな問題を解決することを意図していない。この方法論は、サンプルの1つまたはサブセットに存在する突然変異をより効果的に同定することができうる、非常に密接に関連するサンプルを利用する前提で設計されている。

 

Method

 2つ以上のゲノムを直接比較する1つの方法は、データセットをde novoでアセンブルし、全ゲノムアライメントが可能なMUMmerなどで比較することである。しかし全ゲノムアセンブリは計算コストが高く、かつ誤ったアセンブリを生成する可能性もあり、アライメントによるアプローチよりもさらに大きな問題を引き起こす可能性がある。代わりにDiamundは、ある固定サイズのkについて、すべてのリードで長さkの全シーケンスを数え、次にこれらのk-mersを互いに比較する直接的なアプローチを使用する。アルゴリズムの10の主要なステップを概説する。

 

1、jellyfishを使ったk-mer抽出と頻度カウント。シーケンスしたリードが全て正しければ、(偶発的なマッチがほぼ起こらない)大きめのサイズでk-merをカウントすると、異なるk-merの数はゲノムのサイズ(ここではexomeのサイズである50-60Mb) と等しくなる (ゲノムにリピート領域がない時)。しかし実際は、シーケンスエラーにより異なるk- merの数は爆発的に増える。これにstep2で対処する。

 ↓ 

2、シーケンスエラー由来k-merのフィルタリング。リードに1つのシーケンスエラーが発生すると、本来ゲノムに無いだろうk個の新しいk-mer配列が出来る。そのようなk-merの出現回数は大抵1だが、exome解析ではカバレッジは100近くあり、偶然2以上の出現頻度になるシーケンスエラー由来のk-merは、ゲノム全体では無視できない数になる。著者らのempericalな観測から、頻度10以下のk-merを除くことで、真の情報の損失は限りなくゼロにしつつ、シーケンスエラーはほぼ完全に排除できると判断した。頻度10以下のk-merを捨てる。

  ↓ 

3、ベクターとリファレンスヒトゲノムに出てくるk-merのフィルタリング。ヒトリファレンスゲノムのk-merを除く。これは、リファレンスゲノムに存在するk-mer配列は病気の発生に関与しないと考えられるためである。ベクターデータはUniVec database を使っている。

  ↓ 

4、Proband (子供) と父親、probandと母親それぞれの共通/非共通のk-merの探索。両親のシーケンスデータには見つからないk-merを抽出し、父との比較、母との比較、得られた二つの結果をミックスする。得られたk-merには、probandのexomeに固有のk-merが全て含まれているはずである。

  ↓ 

5、固有のk-merをリードにバックアライメントして、リードに差し戻す。ステップ4の終了時点で情報は1万倍以上濃縮されており、論文のデータでは、数千〜数万のk-merまで減っている(論文 表1)。このk-merを含むリードを抽出し、リードに戻す。MUMmerかmodifyしたkrakenを使いシーケンスリードとのアライメントを行うことで、該当リードを見つけ出す。デフォルトではより高速なkrakenを使う。

  ↓ 

6、再度ベクター配列のフィルタリング。probandsにまだ残っているベクターデータを除く。ステップ3と同じデータベースを使うが、精度の高いvecscreen programを使っている。初期インプットと比べリード数は極端に減っているので、この段階では、コンピュータリソース負荷が大きいプログラムの使用は容易になっている。

 

7、選抜したリードのアセンブリMinimoアセンブラを使い、リードをアセンブルする。ここまでの処理で残ったリードは、変異を含む部位に限定されているはずで、アセンブルしても伸びることはない(変異が集中しない限り)。100-bpのリードならば、アセンブルしても、理屈上、最大199-bpになる。実際はもっと短くなる。

  ↓ 

8、コンティグの選抜。シーケンスエラーにより、同じポジションでもリードにバリエーションが生じ、複数コンティグが生じることがある。bowtieを使ってコンティグをリファレンスゲノムにマッピングする。ゲノムの同じ部位に複数のコンティグがアライメントされるならば、一つのコンティグに統合する。100-bp以下の配列は、100-bpまでリファレンス情報を使い拡張する。(リファレンスゲノムはあくまでコンティグのクラスタリングのためにしか利用していない)。

  ↓ 

9、コンティグの選抜2。両親の両方または片方のexome キャプチャに失敗して、probandに固有の配列に見えている可能性がある。bowtieを使って両親のシーケンスリードを8のコンティグにマッピングして、両親のデータともマッピングされるか確かめる(この時、真のマッピング部位か検証するため、マッピングされたリードは、リファレンスゲノムにもマッピングし、ベストマッピングかどうか確認している)。片親データしかマッピングされなければ、実験のエラーの可能性が捨てきれないので排除する。

またstep8でリファレンスゲノムにアライメントされなかったcontigは、別に保存される(情報がないのでマニュアルで解析する必要がある)。

  ↓ 

10、バリアントのコールとアノテーションsamtoolsのmpileupを使い、リファレンスゲノムと比較して変異している部位を、アミノ酸変化の注釈付きで出力する。

f:id:kazumaxneo:20180410123520j:plain

論文より転載。

 

 

公式サイト

http://ccb.jhu.edu/software/diamund/ 

インストール 

cent OS6に導入した。

 依存

  • Kraken a system for assigning taxonomic labels to short DNA sequences;
  • Samtools package for managing large short read alignment files;
  • Jellyfish a tool for fast, memory-efficient counting of k-mers in DNA;
  • Bowtie2 an ultrafast and memory-efficient tool for aligning sequencing reads;
  • Minimo a de novo assembler based on the AMOS infrastructure;
  • Mummer a system for rapidly aligning entire genomes.
  • vecscreen
  1. Minioは上記リンクからダウンロードして、"./configure"、"sudo make install"する(binにMinimoのバイナリができる)perlのライブラリも必要になるので、cpanm XML::Parserで導入する(DBIStatistics::Descriptiveも同様)。
  2. krakenとjellyfishは以前の説明を参照(リンク)。jellyfishはversion1をインストールすることに注意する。version1ダウンロード方法はkrakenの解説を参照。
  3. 他のツールはbrewで導入できる。
  4. <追記> SAMToolsのコードが旧式のものになっている。旧式のSAMtoolsバイナリをダウンロードして(リンク)、そのSAMtoolsを指定する。 mpileupを使うので、version 0.1.13あたりを使う。perlのコードを現在のSAMtoolsの形式に修正してもよいが、いくつかのサブルーチン内をいじるので手間になる。ソースをダウンロードしてmakeすればsamtoolsができるので、ここにdiamond.cfgのsamtoolsパスを通す。現行のsamtoolsと使い方が少しだが致命的に異なるコマンドがあるので、旧式を/usr/local/bin/samtoolsなどと置き換えてはいけない(sortのコマンドとかmpileupのオプションの違いとか)。
  5. <追記> vecscrrenは公式には記載がないが必要になる。ダウンロードできるURLが見つからなかったので、Githubのvecscreen_plus_taxonomyレポジトリのscritsにあるvecscreenのバイナリをダウンロードした(linux版)(ただしこのvecscreenはバージョンが違うので、Diamondのコードを修正する必要がある)。

    vecscreen_plus_taxonomy/scripts at master · aaschaffer/vecscreen_plus_taxonomy · GitHub

 

本体は公式からダウンロードして解凍する。

tar -xvf diamund.tar
cd Diamund/
perl diamund.pl 

$ perl diamund.pl 

Usage:

 diamund.pl {<patient>.fa|<patient>.fq[.gz|.gz2][,<patient_mates>.fq[.gz|.gz2]]} {<parent1>.fa|<parent1>.fq[.gz|.gz2][,<parent1_mates>.fq[.gz|.gz2]]}  [{<parent2>.fa|<parent2>.fq[.gz|.gz2][,<parent2_mates>.fq[.gz|.gz2]]}] {-e <exon_file>.fa|-x <exon_kmer_file>} [options]

Options:

 -e <exon_file>         fasta file with the sequence of all exons

 -x <exon_kmer_file>    sorted exome kmer file (either the -e or -x options should be present in the input)

 -k|kmer <kmer_size>    kmer size for jellyfish program (default: 27)

 -r|recessive           run program in recessive mode (default: dominant)

 -t <no_of_threads>     number of threads for the jellyfish program

 -v <kmer_vector_file>  file with vector kmers (if available)

 -z <vector_file>       fasta file with vectors

 -f <filter_threshold>  only keep kmers in patient that appear more than the <filter_threshold> times

 -n <read_no>           only keep assemblies that have at least <read_no> assembled reads 

 -d <vecscreen_data>    screen for vector contamination with vescreen using the data from this directory

 -g <[0|1](,[0|1])*>    a string of 0 and 1 separated by commas, and signifying the sex associated with each data sample respectively; 1 is male, and 0 is female  

 -l <read_len>          read length (default:100)

 -o <output_file>       output file with variants called (default is STDOUT)

 -s <sequence_depth>    observed depth of coverage in non-patients in order to call variant (default:10)

 -y <kmer_no>           assume kmers occuring this many times in normals are errors (default: 0)

 -c                     don\'t clean-up the files                  

 -m                     use mummer instead of kraken to align kmers to reads (default: no)

Usage examples:

 diamund.pl patient.fa father.fa mother.fa -e hg19_exome.fa -g 1,0,1

or

 diamund.pl patient.1.fq.gz,patient.2.fq.gz normal.1.fg.gz,normal.2.fg.gz -x 27mers-from-exome-sort

Error: at least two input fasta file need to be specified

 

 公式からhuman hg19のFASTAとbowtie indexをダウンロードする。

wget ftp://ftp.ccb.jhu.edu/pub/software/diamund/data/hg19.tar.gz
mkdir human_ref_index
tar -zxvf hg19.tar.gz -C human_ref_index/ #human_ref_index/に解凍

Vectorとexomeの計算済みk-merデータベースも公式に準備されている。ダウンロードする。

wget ftp://ftp.ccb.jhu.edu/pub/software/diamund/data/27mers_exome_vect.tar.gz
mkdir k-mer_databae
tar -zxvf 27mers_exome_vect.tar.gz -C k-mer_databae/

ベクターVeccscreenのデータベース配列へのリンクも公式に載っているが、NCBI ftpサーバーのURLが変更されている(The UniVec Database)。以下のURLでダウンロードできる。

wget ftp://ftp.ncbi.nlm.nih.gov/pub/UniVec/UniVec

UniVecはFASTAフォーマットのデータ。

$ head UniVec

>gnl|uv|X66730.1:1-2687-49 B.bronchiseptica plasmid pBBR1 genes for mobilization and replication

CTCGGGCCGTCTCTTGGGCTTGATCGGCCTTCTTGCGCATCTCACGCGCTCCTGCGGCGGCCTGTAGGGC

AGGCTCATACCCCTGCCGAACCGCTTTTGTCAGCCGGTCGGCCACGGCTTCCGGCGTCTCAACGCGCTTT

GAGATTCCCAGCTTTTCGGCCAATCCCTGCGGTGCATAGGCGCGTGGCTCGACCGCTTGCGGGCTGATGG

TGACGTGGCCCACTGGTGGCCGCTCCAGGGCCTCGTAGAACGCCTGAATGCGCGTGTGACGTGCCTTGCT

GCCCTCGATGCCCCGTTGCAGCCCTAGATCGGCCACAGCGGCCGCAAACGTGGTCTGGTCGCGGGTCATC

TGCGCTTTGTTGCCGATGAACTCCTTGGCCGACAGCCTGCCGTCCTGCGTCAGCGGCACCACGAACGCGG

TCATGTGCGGGCTGGTTTCGTCACGGTGGATGCTGGCCGTCACGATGCGATCCGCCCCGTACTTGTCCGC

CAGCCACTTGTGCGCCTTCTCGAAGAACGCCGCCTGCTGTTCTTGGCTGGCCGACTTCCACCATTCCGGG

blastのデータベースを作る。

makdir database
makeblastdb -dbtype nucl -in UniVec -out database/UniVec
ls database/

$ ls database/

UniVec.nhr  UniVec.nin  UniVec.nsq

 

 diamond.cfgを開き、各ツールのパスを記載しておく。さきほどダウンロードしたデータベースのパスと、コマンドのパスから(whichで調べる)、以下のように変更した。

$ cat diamund.cfg 

mummer  ~/.linuxbrew/bin/mummer

vecscreen_data  /home/uesaka/Diamund/database/UniVec

kraken  /home/disk2/uesaka/kraken

human_index /home/uesaka/Diamund/human_ref_index/hg19

human_fasta /home/uesaka/Diamund/hg19.fa

jellyfish /home/disk2/uesaka/kraken/jellyfish-1.1.11/bin/jellyfish

minimo  ~/amos-3.1.0/bin/Minimo

bowtie  ~/.linuxbrew/bin/bowtie

bowtie2  ~/.linuxbrew/bin/bowtie2

samtools /home/uesaka/Diamund/samtools-0.1.13/samtools

bcftools  /usr/local/bin/a5_miseq_linux_20160825/bin/bcftools

vecscreen /user/local/bin/vecscreen

 bowtieのほか、はbowtie2の行を追加する。 

 

vecscrren2を使ったので、オプションを変更する必要があった。結果が変わる可能性を覚悟で988行目の

my $out= `echo $seq | $vecscreen -d $vecscreen_data -f3`;

my $out= `echo $seq | $vecscreen -db $vecscreen_data`;

に変更した。

 

 

 

 

ラン

 fammily triosの解析。3人のシーケンスデータを入力する。

diamund.pl patient.fa father.fa mother.fa -e hg19_exome.fa -g 1,0,1 
  • -e <exon_file>     fasta file with the sequence of all exons 
  • -g <[0|1](,[0|1])*>   a string of 0 and 1 separated by commas, and signifying the sex associated with each data sample respectively; 1 is male, and 0 is female

正常組織と病原組織の解析。fastqはgz圧縮形式にも対応している。

diamund.pl patient.1.fq.gz,patient.2.fq.gz normal.1.fg.gz,normal.2.fg.gz -x 27mers-from-exome-sort
  •  -x <exon_kmer_file>   sorted exome kmer file (either the -e or -x options should be present in the input)

jellyfishはメモリが足りないとsegmentation errorを起こす。WESでも複数走らせる時は数百GBいるかもしれない。注意する。

 

 

テストデータとして、 Manuel Corpasさんが公開してくれている親子トリオのexomeデータを使ってみる。

https://cambridgeprecisionmedicine.com/2011/07/12/getting-my-genome-sequencing-done-part-i/

https://cambridgeprecisionmedicine.com/2013/01/21/corpas-family-exome-data-available-for-public-download/amp/

slideshareからダウンロードし解凍する。

unzip 106340.zip

dad_1.clean.fq.gz、dad_2.clean.fq.gz、mom_1.clean.fq.gz、mom_2.clean.fq.gz、daughter_1.clean.fq.gz、daughter_2.clean.fq.gzができる。

父、母、娘のデータを指定してランする。

perl diamund.pl dad_1.clean.fq.gz,dad_2.clean.fq.gz mom_1.clean.fq.gz,mom_2.clean.fq.gz daughter_1.clean.fq.gz,daughter_2.clean.fq.gz -s 10 -t 40 -x k-mer_databae/27mers-exome

 

krakenを使ってk-merからリードに戻すstepでエラーになる。ツールのバージョン違いの気が下が、手っ取り早く-mをつけてkrakenからmummerに切り替えて対応した。

bcfしか出力しなかったので、bcftoolsを使いvcfに変更した。

bcftools view 2-to-1.bcf > out.vcf

個人データの結果は載せないが、出力は確かにかなり限定された。

 

 

 

感想

思いっきりempericalなアルゴリズムwiki)に頼った手法なので、条件を満たせば良い結果を出してきます。数学的な拠り所が少ないのでインフォマティクスが専門の方には興味が持たれないツールかもしれませんが、wetの研究者にとっては良いツールではないでしょうか?

 

ただし、このツールはどうやら堅牢性があまり考慮されておらず(本体はperlスクリプト)、ツールのバージョンが少しでもずれるとエラーを起こし、そこでジョブが止まります。マニュアル不足で3rd partyのツールのインストールに手間取ったのも影響してますが、繰り返しテストしていたら1週間以上かかりました。ランニング時は > log.txt 2>&1とかをつけて、標準出力のエラーログを取りながら、どのstepでつまづいたか1つ1つ確認して解決していくのが良いでしょう。またテスト時はchr1の一部などにマッピングして、マッピングされたfastqだけ回収し、それで進めてください。速度が考慮されておらず、テストにもかなり時間がかかります。

 

現在開発中となってますが、RUFUSも発想は同様のツールのようです(ただしユニークk-merを抽出してそれをそのままアセンブリする)。

GitHub - jandrewrfarrell/RUFUS: RUFUS k-mer based genomic variant detection

 

 

引用

DIAMUND: Direct Comparison of Genomes to Detect Mutations

Steven L Salzberg, Mihaela Pertea, Jill A Fahrner, and Nara Sobreira

Hum Mutat. 2014 Mar; 35(3): 283–288.

 

RNA seqのシミュレータ polyester

 

 RNA-seq実験は遺伝子発現を研究する手段としてますます普及が進んでいる。RNA-seqデータ(Oshlack et al、2010)の発現解析のための様々な統計的手法がある。 RNA-seqの統計的方法論の開発者は、ツールが正しく機能しているかどうかをテストする必要がある。真の遺伝子発現レベルおよび集団間の発現差異は通常不明であり、スパイクイン実験は時間と経費の両方でコストがかかるため、実際のデータセットで精度試験を行うことはできない。

 代わりに、研究者は計算シミュレーションを使用して、既知の信号とノイズ構造を持つデータセットを作成することがよくある。典型的には、発現変動分析ツールを評価するために使用されるシミュレートの発現は、一般的な発現変動分析ツール(AndersおよびHuber、2010; Robinson et al、2010)において使用される統計モデルから遺伝子カウントとして生成される。しかし、これらのシミュレートされたシナリオは、シーケンスリードまたはリードカウントのようなRNA-seqデータ分析の上流工程中に生じる発現の変動性を説明していない。polyesterは、RNA-seqのリードをシミュレートするための新しいRパッケージである。polyesterの主な利点は、ユーザーが、遺伝子またはアイソフォームの特定の発現変動シグナルを用いてシークエンシングリードをシミュレートできることである。これにより、ユーザーはRNA-seqパイプラインの複数の観点で変動の原因を調べることができる。

 シーケンシングリードを生成する既存のRNA-seqシミュレータは、実際のレプリケート実験および指定された発現変動シグナルを用いた実験のために設計されていない。たとえば、RSEM(Li and Dewey、2011)に同梱されているrsem-simulate-readsユーティリティでは、実際のシークエンシングリードをアライメントさせてモデルを開発するため、シミュレーションに時間がかかる。FluxSimulator(Griebel et al、2012)もRNA-Seqの有効性を評価するためのBenchmarkerであるBEERS(Grant et al、2011)も、発現変動を誘導する組み込みメカニズムを持っていない。これらのシミュレーターは、レプリケート間の生物的変動や特定の転写物の正確な発現レベルを特定するための方法も提供しない。 TuxSimは、微分式(Trapnell et al、2013)を用いてRNA-seqデータセットをシミュレートするために使用されているが、公に利用可能ではない。

 著者らは、レプリケートおよび明確に定義された発現変動でのRNA-seqのリードをシミュレートするためにpolyesterを作製した。ユーザーは、少数の遺伝子または単一の染色体から小さな実験を簡単にシミュレートできる。これによって、計算集約的なステップを実行する必要がある場合に、シミュレーションの計算時間を短縮することができる。

 

 

インストール

source("http://bioconductor.org/biocLite.R") 
biocLite("polyester")

 

ラン

シミュレートにはトランスクリプトFASTAファイルが必要である。chr20のtranscriptsが用意されている(Users/user_name/Library/R/3.4/library/polyester/extdata/)。これを使ってみる。

 

ライブラリの読み込み。

library(polyester) 
library(Biostrings)

 

transcriptsのFASTAファイルの読み込み。

fasta_file = system.file('extdata', 'chr22.fa', package='polyester') #パッケーージに入っているファイルを探すsystem.file関数
fasta = readDNAStringSet(fasta_file) #BiostringsのreadDNAStringSet

> head(fasta,3) 

> head(fasta,3)

  A DNAStringSet instance of length 3

    width seq                                               names               

[1]   725 GCCACGTGAAGGATGTGTTTGCT...TTAAATTTGATTGAACAATTGTA gi|424037187|ref|...

[2]   746 GCCACGTGAAGGATGTGTTTGCT...TAAAAAAAAAAAAAAAAAAAAAA gi|424037186|ref|...

[3]  2037 GGTAGACGCGATCTGCTGGCTAC...GAACATGATGGGGGGCATGCGTC gi|209977002|ref|...

 

テストデータは900以上のtranscriptsがあるが、先頭20配列だけ使う。 

small_fasta = fasta[1:20]

#ファイル名chr22_small.faで書き出す。このあと使う。
writeXStringSet(small_fasta, 'chr22_small.fa')

 

全transcripts x20のFPKMとする(20xカバレッジ = transcriptsサイズ / リード長 x 20)

readspertx = round(20 * width(small_fasta) / 100)

 

matrix関数で行列を20行作成させる。これが次のステップのsimulate_experiment関数に与える各transcriptsのfold_changesの値になる。

fold_changes = matrix(c(4,4,rep(1,18),1,1,4,4,rep(1,16)), nrow=20) #rep(1,18)は1を18回繰り返し
head(fold_changes,20)

> head(fold_changes,20)

      [,1] [,2]

 [1,]    4    1

 [2,]    4    1

 [3,]    1    4

 [4,]    1    4

 [5,]    1    1

 [6,]    1    1

 [7,]    1    1

 [8,]    1    1

 [9,]    1    1

[10,]    1    1

[11,]    1    1

[12,]    1    1

[13,]    1    1

[14,]    1    1

[15,]    1    1

[16,]    1    1

[17,]    1    1

[18,]    1    1

[19,]    1    1

[20,]    1    1

20行だと、1列目は4,4,rep(1,18)の20文字。2列目は1,1,4,4,rep(1,16)の20文字。

 

負の二項分布でシミュレートするsimulate_experiment関数でリードを発生させる(リンク)。サンプル1もサンプル2も3のレプリケートデータを発生させるc(3,3)。

simulate_experiment('chr22_small.fa', reads_per_transcript=readspertx, num_reps=c(3,3), fold_changes=fold_changes, outdir='simulated_reads', paired=FALSE)
  • reads_per_transcript  baseline mean number of reads to simulate from each transcript. Can be an integer, in which case this many reads are simulated from each transcript, or an integer vector whose length matches the number of transcripts in fasta. (Default 300)
  • fold_changes  Vector of multiplicative fold changes between groups, one entry per transcript in fasta. A fold change > 1 means the transcript is overexpressed in the first num_reps (or num_reps[1]) samples. Fold change < 1 means transcript is overexpressed in the last num_reps (or num_reps[2]) samples. The change is in the mean number of reads generated from the transcript, between groups.
  • num_reps    How many biological replicates should be in each group? If num_reps is a single integer, num_reps replicates will be simulated in each group. Otherwise, num_reps can be a length-2 vector, where num_reps[1] and num_reps[2] replicates will be simulated in each of the two groups.
  • paired   If TRUE, paired-end reads are simulated; else single-end reads are simulated.

出力フォルダ。

f:id:kazumaxneo:20180408205902j:plain

ペアエンドならsample_01_1.fasta、sample_01_2.fastaなどの名前になる。

 

リード数を直接指定する例は、Bioconductorのオンラインマニュアルを確認してください。オンラインマニュアルには、実際の発現データの行列(FPKMまたはRPKM)およびアノテーションを取り込んで、与えられたデータと同様の発現レベルおよび発現変動をシミュレートする関数simulate_experiment_empiricalについても説明されています。

 

引用

Polyester: simulating RNA-seq datasets with differential transcript expression

Alyssa C. Frazee, Andrew E. Jaffe, Ben Langmead, and Jeffrey T. Leek.

Bioinformatics. 2015 Sep 1; 31(17): 2778–2784. 

 

低複雑度領域由来のリードを除去する RepeatSoaker

 

 次世代シークエンシング(NGS)技術は、主に、DNA / RNAサンプルからの数百万回のリードの超並列シーケンシングに基づいており、リード長は増加している[論文より ref.1,2]。 NGSのコストは急速に低下し、その結果、転写を研究するためにマイクロアレイの使用からRNA-seqデータへの移行が比較的急速に進んでいる。このNGSへの依存度の高まりは、データの品質やデータの解釈に影響を及ぼす可能性のある分析ステップの検討を必要とする。

 さまざまなタイプのNGS実験やライブラリー作成プロトコル下流の処理ステップを決定しているが、ライブラリを構築するために使用したアダプター配列を除去し、低クオリティな塩基を除去する[ref.5]。これに続いて、PCRによるライブラリー増幅の際に起こり得る、重複したリードを除去する[ref.6]。このステップの背後にある論理的根拠は、このような重複したリードは、DNA-seqデータ[ref.7]における変異検出、RNA-seqデータにおける遺伝子発現[ref.8]、およびChIP-seqデータ[10]において遺伝子の定量化など生物学的シグナルの真のレベルに関する誤った結論を起こすということである[ref.9]。この問題にどのように最善を尽くすかについて、2つの考え方が現場で浮上している。第1は、潜在的なバイアスだと仮定して、データセットからすべての重複またはlow complexity (以後、低複雑)なリードを除去する[ref.7,8,10,11,12]。二つ目は、これらの複製が真の固有の観測であり、それらの除去がそれ自身のバイアスを導入すると信じている[ref.7]。重複リードの影響はDNAシークエンシングで研究されているが、ChIP-seq実験におけるモチーフ検出およびRNA-seq実験における遺伝子発現による重複がどのように遺伝子シグナルに影響を与えるかという問題は未解決のままである。

 リファレンスゲノムにおける低複雑性[ref.14,15]および反復[ref.16]の存在は、それらが生物学的研究の結論にどのように影響を及ぼすかにあまり注意を払わなかった。そのような領域は、それらから発生するリードが複数の場所にマッピングされ、位置合わせが複雑になる。この問題は、真核生物ゲノムがリピート領域において非常に豊富であり得るため、小さくはない。例えば、ヒトゲノムには〜47%のリピートが含まれていると推定している研究者もいる[ref.19]。最近のENCODEプロジェクトの結果は、ゲノムが広く転写されていることを示唆しているが[ref.20、21]、低複雑性およびリピート領域に由来するRNA分子は、ゲノム全体にわたって無差別にアライメントされる可能性がある(一部略)。

 この研究では、RNA-seqおよびChIP-seq実験において、アダプタートリミング、重複リード除去、および低複雑性領域に重複するリードを排除する効果を系統的に調査した。各処理段階で、研究者は、RNA-seqデータから検出された遺伝子の経路および遺伝子のオントロジー濃縮分析およびChIP-seqデータから検出されたピークにおける転写因子結合部位濃縮分析を実施することによって、生物学的シグナルの強度の間接的測定を使用した。ここでの著者らの理論的根拠は、処理ステップがより有意な濃縮p値をもたらす場合、その処理ステップは生物学的シグナルに積極的に影響を及ぼす可能性があるということである。

 低複雑度領域をオーバーラップするリードの除去は、それほど注意を払われなかった。例えば、このステップは、PRINSEQツール[ref.25]の品質管理ステップのセットに含まれている。低複雑度合いの影響を個別にテストできるようにするために、簡単なポストアライメントフィルタリングツール、RepeatSoakerを使用して、複雑さの低い領域のゲノム座標を含むユーザー提供のテンプレートファイルと重複する読み取りを除外する。アライナーに依存しないように設計されたRepeatSoakerは、アライメントされたデータをBAM形式で処理し、複雑さの低いリードを取り除き、クリーンなBAMファイルとフィルタリング統計を出力する。RepeatSoakerは、NGSデータからアライメントアーチファクトを除去する簡単な方法で、転写物発現の定量化における誤った陽性の可能性を排除するように設計されている。低複雑な領域を含むリードがChIP-seqのようなバイアスを導入する可能性がある他のシーケンシング技術に拡張可能であり、著者らはRepeatSoakerが再現性のあるNGSパイプラインをよりよく構造化するのを助ける標準的なステップになると考えている。

 著者らは、アダプタートリミング、重複除去、およびRNA-seqおよびChIP-seq実験に対するRepeatSoakerを用いた複雑なリードのフィルタリングを適用し、各ステップが下流の濃縮分析の結果にどのように影響するかを調べた。結果は、アダプタートリミングが、RNA-seqデータにおける遺伝子オントロジーおよび経路濃縮分析の重要性を高め、ChIP-seqデータにおけるモチーフ検出を強化することを示していた。リード数を減らしたにもかかわらず、重複除去ステップは、特にアダプタートリミングと組み合わせた場合、生物学的シグナルの重要性を高めるのにさらに役立った。 RepeatSoakerを使用して複雑さの低いリードをフィルタリングすることは、リード総数にはあまり影響しないが、このステップは生物学的信号の検出に正の効果をもたらした。この研究では、アダプタートリミングと重複除去は、RNA-seqとChIP-seqデータ内でより強い生物学的シグナルを検出する上で重要なステップであることを示唆していた。

 

このRepeatSoaker自身に低複雑度の領域を探す機能はない。RepeatMaskerで低複雑度の領域を探し、RepeatSoakerにbedで与えると、入力したbamからその領域にアライメントされたリードが除外される。

 

インストール

依存

  • bedtools

Github

https://github.com/mdozmorov/RepeatSoaker

git clone https://github.com/mdozmorov/RepeatSoaker.git
make mappability
make clean

> ./repeat-soaker

$ ./repeat-soaker 

USAGE: repeat-soaker <options> -r repeat_regions.bed in.bam

 

Options:

    -o <out.bam> : specify an output path for the bam

         default: in.repeatSoaker.bam

    -h : show this help

    -p : percent overlap (a float between 0 and 1)

         default: 0.75

 

 ラン

ダウンロードしたフォルダには、人の low complexilityな領域の1例が含まれる。新規に探索するには、RepeatMasker(-noint/-intつき)で low complexilityの領域を検出し、BEDOPS rmsk2bedでRepeatMasker出力をbedに変換する。

 

除去したい領域のbedを指定する。 入力のbamはcoordinate ソートされている必要がある。samtools sortコマンドを使う。

samtools sort -@ 8 input.bam > coordsorted.bam

#またはpicard-toolsを使う。
picard SortSam INPUT=input.bam OUTPUT=coordsorted.bam SORT_ORDER=coordinate

 

除去したいlow complexility領域のbedを指定してランする。

repeat-soaker -r rmsk.bed -o output.bam coordsorted.bam

 

 

引用

Detrimental effects of duplicate reads and low complexity regions on RNA- and ChIP-seq data.

Dozmorov MG, Adrianto I, Giles CB, Glass E, Glenn SB, Montgomery C, Sivils KL, Olson LE, Iwayama T, Freeman WM, Lessard CJ, Wren JD.

BMC Bioinformatics. 2015;16 Suppl 13:S10.

 

トリミングツール fqtrim

 

 fqtrimは、アダプター、polyA tail、未知塩基(Ns)および低クオリティな3 '領域をトリミングできる多目的トリミングツール。アダプターとポリA配列の不正確なマッチングにも対応している。 このユーティリティは、複雑さの低い配列(ダスト)のフィルターを適用したり、重複リードをカウントすることもできる。例えばこの論文(リンク)では、low complexilityなONTリードの除去に使用されている。

 

公式サイト

https://ccb.jhu.edu/software/fqtrim/

 

インストール

Github

https://github.com/gpertea/fqtrim

公式よりダウンロードするか、Githubからcloneしてビルドする。

git clone https://github.com/gpertea/fqtrim.git
cd fqtrim/
make release

> ./fqtrim

$ ./fqtrim 

fqtrim v0.9.7. Usage:

fqtrim [{-5 <5adapter> -3 <3adapter>|-f <adapters_file>}] [-a <min_match>]\

   [-R] [-q <minq> [-t <trim_max_len>]] [-p <numcpus>] [-P {64|33}] \

   [-m <max_percN>] [--ntrimdist=<max_Ntrim_dist>] [-l <minlen>] [-C]\

   [-o <outsuffix> [--outdir <outdir>]] [-D][-Q][-O] [-n <rename_prefix>]\

   [-r <trim_report.txt>] [-y <min_poly>] [-A|-B] <input.fq>[,<input_mates.fq>\

 

 Trim low quality bases at the 3' end and can trim adapter sequence(s), filter

 for low complexity and collapse duplicate reads.

 If read pairs should be trimmed and kept together (i.e. never discarding

 only one read in a pair), the two file names should be given delimited by a comma

 or a colon character.

 

Options:

-n  rename the reads using the <prefix> followed by a read counter;

    if -C option was also provided, the suffix "_x<N>" is appended

    (where <N> is the read duplication count)

-o  write the trimmed/filtered reads to file(s) named <input>.<outsuffix>

    which will be created in the current (working) directory (unless --outdir

    is used); this suffix should include the file extension; if this extension

    is .gz, .gzip or .bz2 then the output will be compressed accordingly.

    NOTE: if the input file is '-' (stdin) then this is the full name of the

    output file, not just the suffix.

--outdir for -o option, write the output file(s) to <outdir> directory instead

-f  file with adapter sequences to trim, each line having this format:

    [<5_adapter_sequence>][ <3_adapter_sequence>]

-5  trim the given adapter or primer sequence at the 5' end of each read

    (e.g. -5 CGACAGGTTCAGAGTTCTACAGTCCGACGATC)

-3  trim the given adapter sequence at the 3' end of each read

    (e.g. -3 TCGTATGCCGTCTTCTGCTTG)

-A  disable polyA/T trimming (enabled by default)

-B  trim polyA/T at both ends (default: only poly-A at 3' end, poly-T at 5')

-O  output only reads affected by trimming (discard clean reads!)

-y  minimum length of poly-A/T run to remove (6)

-q  trim read ends where the quality value drops below <minq>

-w  for -q, sliding window size for calculating avg. quality (default 6)

-t  for -q, limit maximum trimming at either end to <trim_max_len>

-m  maximum percentage of Ns allowed in a read after trimming (default 5)

-l  minimum read length after trimming (if the remaining sequence is shorter

    than this, the read will be discarded (trashed)(default: 16)

-r  write a "trimming report" file listing the affected reads with a list

    of trimming operations

-s1/-s2:  for paired reads, one of the reads (1 or 2) is not being processed

    (no attempt to trim it) but the pair is discarded if the other read is

    trashed by the trimming process

--aidx option can only be given with -r and -f options and it makes all the 

    vector/adapter trimming operations encoded as a,b,c,.. instead of V,

    corresponding to the order of adapter sequences in the -f file

-T  write the number of bases trimmed at 5' and 3' ends after the read names

    in the FASTA/FASTQ output file(s)

-D  pass reads through a low-complexity (dust) filter and discard any read

    that has over 501506411212f its length masked as low complexity

--dmask option is the same with -D but fqtrim will actually mask the low 

    complexity regions with Ns in the output sequence

-C  collapse duplicate reads and append a _x<N>count suffix to the read

    name (where <N> is the duplication count)

-p  use <numcpus> CPUs (threads) on the local machine

-P  input is phred64/phred33 (use -P64 or -P33)

-Q  convert quality values to the other Phred qv type

-M  disable read name consistency check for paired reads

-V  show verbose trimming summary

Advanced adapter/primer match options (for -f or -5 , -3 options):

  -a      minimum length of exact suffix-prefix match with adapter sequence that

          can be trimmed at either end of the read (default: 6)

  --pid5  minimum percent identity for adapter match at 5' end (default 96.0)

  --pid3  minimum percent identity for adapter match at 3' end (default 94.0)

  --mism  mismatch penalty for scoring the adapter alignment (default 3)

  --match match reward for scoring the adapter alignment (default 1)

  -R      also look for terminal alignments with the reverse complement

          of the adapter sequence(s)

 

ラン

クオリティ30以下の部位をトリミングし、25bp以下になったリードは除去する。

fqtrim -q 30 -l 25 -o filtered.fq input.fq
  • -P    input is phred64/phred33 (use -P64 or -P33)
  • -q    trim read ends where the quality value drops below <minq>
  • -l     minimum read length after trimming (if the remaining sequence is shorter than this, the read will be discarded (trashed)(default: 16)
  • -o    write the trimmed/filtered reads to file(s) named <input>.<outsuffix> which will be created in the current (working) directory (unless --outdir is used); this suffix should include the file extension; if this extension is .gz, .gzip or .bz2 then the output will be compressed accordingly. NOTE: if the input file is '-' (stdin) then this is the full name of the output file, not just the suffix.
  • -p    use <numcpus> CPUs (threads) on the local machine
  • -P    input is phred64/phred33 (use -P64 or -P33)

input. filtered.fqが出力される。

 

 low complexilityフィルターをかける。低複雑度の配列が除去される。

fqtrim -D -o filtered.fq input.fq
  •  -D   pass reads through a low-complexity (dust) filter and discard any read that has over 501577251212f its length masked as low complexity
  • --dmask   option is the same with -D but fqtrim will actually mask the low 

 

引用

https://ccb.jhu.edu/software/fqtrim/

 

https://zenodo.org/record/1185412#.WsmDTSPAORc

 

NGSデータから素早くバクテリアの分析を行う MICRA

 

 

 ハイスループットシーケンシング(HTS)技術は多くの微生物学的問題に対処するための費用対効果の高い便利なアプローチとして浮上し、この分野を大きく変えている。完全なゲノム情報にアクセスすることは、微生物学における基礎研究に革命をもたらし、例えば、新規な抗菌化合物およびワクチンの開発を可能にしている[論文より ref.1]。 HTSはまた、ゲノムベースの診断[ref.2]のための新しい方法を提供し、したがって感染症[ref.3]およびアウトブレイクのフォローアップ[ref.4]の治療のための臨床応用のための新しいアプローチを可能にする。微生物ゲノムを調査するためのスループットの点で十分なパフォーマンスを備えた安価なプラットフォームである第2世代のベンチトップシーケンサーの発売は、微生物学におけるこの革命に貢献しており、第3世代のシーケンサーはこれを加速する。食品医薬品局(FDA)は、感染症診断と抗菌剤耐性および病原性マーカー検出のためのHTS検査の検証のためのコンセプトを開発しており、HTSベース診断のマイルストーンとなっている。しかし、実験室で実施されるためには、技術的側面とバイオインフォマティクス分析の両方に対して、標準化、最適化、および検証アッセイを実施しなければならない。

 膨大な数のツールが利用可能であるが、HTSデータを分析することは、基礎的および臨床的研究にとっては困難なことである。技術とアルゴリズムの絶え間ない進化は、HTS分析のための標準化された手続きの欠如をもたらす[ref.6]。さらに、バイオインフォマティクスがサポートされていないと、完全な分析のためにいくつかのツールを連続して使用する必要がある。たった数個の完全なパイプラインが存在し、それらの大部分は病原体同定のためだけである。 「識別」という用語は、「特徴付け」という用語と区別されなければならない。同定の目的は菌株を決定することであるが、微生物ゲノムの特徴付けは、遺伝子移入または特定の変異を含む遺伝的特徴を強調する深い研究である。微生物ゲノムの特徴付けのための使いやすいインターフェースと判読可能な結果を​​提供する迅速で自動化されたツールは、HTSデータを時間的に完全に活用する上で非常に重要になる。

 この論文では、リードの分析による微生物同定と特性評価のためのMICRAについて説明する。 MICRAの独創性は、デノボアセンブリおよび配列アノテーションの普及した手順に従うのではなく、効率的なリードのマッピング方法と登録された微生物ゲノムを利用する方法にある。シーケンシング深度のバラツキ(バクテリアがプラスミドを含む場合など)は、デノボアセンブリプロセスを困難にし、実行時間を大幅に増加させる可能性がある。 MICRAでは2種類のマッパー、すなわち高速なマッパーと堅牢なマッパーが結合され、効果的かつデータの損失を回避する原核生物に特化した新しいバリエーションコーラーが開発された[ref.7]。 MICRAは、臨床医と生物学者の両方にとって使いやすいように開発されたWebインターフェイスとして自由に利用できる[ref.8]。パラメータは、ユーザーに完全な自動分析を提供するように最適化されており、要求されるのはリード情報のみになる。より専門的な作業のために、MICRAは、多くの設定パラメータへのアクセスを提供し、カスタマイズ可能にする。出力は直接読み取って使用可能である。 MICRAはいくつかの研究で評価された。すべての結果は、MICRAが迅速であり(ほとんどの場合約10分)、基本的な微生物学ならびにcirculating and emerging strainsの研究のための新しい洞察を提供することを示している。

 

fastqをアップデートすると、自動でベストマッチするリファレンスを決定し、バリアントを検出する。いくつかの処理でアンマップのままのリードはアセンブルされ、病原性細菌にフォーカスしたデータベースPATRICに対して相同性検索される。また、抗生物質耐性遺伝子の検索も行われる。全ての処理は数分で終わる。高速なマッパーと精度の高いマッパーを使い分けることで、迅速かつ精度の高い分析が可能になっている。

 

f:id:kazumaxneo:20180407194227j:plain

論文より転載。MICRAのパイプライン。

 

シーケンスデータをアップロードして解析できるwebサーバーが準備されている。

http://www.pegase-biosciences.com/MICRA/micra.php

f:id:kazumaxneo:20180407152922j:plain

ユーザーガイド

http://www.pegase-biosciences.com/MICRA/help.php

 

メールアドレスを記載し、fastqファイルを指定する。gzとzip圧縮にも対応している。ファイルサイズは最大1GB。illuminaのほかion torrentのシーケンスデータにも対応している。

f:id:kazumaxneo:20180407185326j:plain

 

advancedにチェックをつける。チェックをいれると、処理の有無や、パラメータ設定などをユーザーが指定できる。

f:id:kazumaxneo:20180407185415j:plain

 

解析前にQCとクオリティトリミングも行える。

f:id:kazumaxneo:20180407185455j:plain

 

デフォルトではMICRAがベストマッチのリファレンスを検索し、それにマッピングを行うが(詳細はマニュアル参照)、ユーザーがリファレンスを提供することもできる。

f:id:kazumaxneo:20180407190224j:plain

 

General options

f:id:kazumaxneo:20180407190307j:plain

 

Performs Antibiotic moduleにチェックを入れると、抗生物質耐性遺伝子を検索できる。

f:id:kazumaxneo:20180407190429j:plain

Jobをスタートする。早ければ数分で解析終了のメールが届く。

 

解凍したフォルダ。複数フォルダがあるがsummary.htmlがまとめとなる。

f:id:kazumaxneo:20180407192722j:plain

summary.html

 

f:id:kazumaxneo:20180407193208j:plain

Image for the comparative genomics (SVG): 使用されたリファレンスと比べて違いがある遺伝子は、色が変化する。

f:id:kazumaxneo:20180407193351j:plain

 

 

 

そのほかMICRA: Comparison moduleでは、MICRAの出力ファイルを使い複数のゲノムを比較してベン図などを描画できる。

MICRA: comparison module

 

論文で使用されたテストデータもダウンロードできる。

http://www.pegase-biosciences.com/MICRA/data.html

このリンクから、Clinical case: the 2011 German outbreak caused by Escherichia coli O104:H4、の解析結果をダウンロードした。DrugBank matchとARDB matchの抗生物質耐性遺伝子検出結果を開いてみる。

f:id:kazumaxneo:20180407184454j:plain

複数の抗生物質耐性遺伝子が検出された。

 

引用

MICRA: an automatic pipeline for fast characterization of microbial genomes from high-throughput sequencing data

Ségolène Caboche, Gaël Even, Alexandre Loywick, Christophe Audebert and David Hot

Genome Biology 2017 18:233