WIGファイルの圧縮と解凍を行う smallWig

　NGSのシークエンシング技術の発達により、DNA / RNAのシーケンスと発現解析のコストが劇的に減少した。 RNA-seqは、様々な種および生物、ならびに異なる器官および細胞集団の全トランスクリプトーム情報を提供する、重要かつ安価な技術になった。RNA-seq実験によって生成された膨大な量のデータは、データ記憶コストおよび通信帯域幅の要求を大幅に増加させる。 bigWigやcWigなどのRNA-seqデータ用の現在の圧縮ツールは、汎用の圧縮器（gzip）を使用しているが、最適な圧縮方式については改善の余地がある。著者らはsmall-wigというWIGデータ用のロスレス圧縮方式を新たに開発した。

　RNA-seq WIGファイルに対してsmallWig圧縮アルゴリズムのさまざまなバリエーションをテストした。その結果、smallWigは平均して18倍の圧縮率向上、最大2.5倍の圧縮時間の改善、bigWigと比較した場合の1.5倍の減圧時間の改善を示している。テストされたファイルでは、アルゴリズムのメモリ使用量は決して90 KBを超えなかった。考案された形式は、視覚化、要約統計解析、圧縮ファイルからの高速照会を可能にするランダムアクセス機能を備えており、bigWigと比較して桁違いの改善をもたらし、圧縮率がcWigによって生成されるもののほんの一部にすぎないことを示している。

インストール

cent OS6に導入した。

Github

git clone https://github.com/brushyy/smallWig.git
cd smallWig/
make 

> ./smallwig2wig

$ ./smallwig2wig 

To use:

 smallwig2wig [InputFile] [OutputFile] 

options: 

-s [ChrmName (e.g. chr1)] [Query Start (integer)] [Query End (integer)] 

subsequence query 

-p [total number of processes]

parallel realization, only availabe if -s is disabled and -r is enabled in encoding

> ./wig2smallwig 

$ ./wig2smallwig 

To use:

wig2smallwig InputFile OutputFile 

options: 

-m [N=0..9, uses 3 + 3*2^N MB memory, decompress should use same N] 

context mixing 

-r [B, encode block size from 8 to 32]

random access and encode by blocks of size 2^B

-p [total number of processes] 

parallel realization, only available if -r is enabled and -m is disabled

ラン

bamからwiggle形式への変換（紹介）(*1 )。

bam2wig input.bam

出力ディレクトリが作られ、中に input.wigができる。 

4スレッドパラレルに使い、smallwigに変換。

wig2smallwig input.wig out.swig -p 4

•  -m [N=0..9, uses 3 + 3*2^N MB memory, decompress should use same N] 

全て解凍する。

smallwig2wig in.swig out.wig

chr1の300-500だけ解凍する。

smallwig2wig in.swig out.wig -s chr1 300 500

引用

smallWig: parallel compression of RNA-seq WIG files.

Wang Z, Weissman T, Milenkovic O.

Bioinformatics. 2016 Jan 15;32(2):173-80.

*1

似たコマンドが他にも存在するので注意してください。

50近いバクテリアの１万以上の機能未知遺伝子欠損の影響をまとめた Fitness Browser

注意: タイトルには 機能未知遺伝子だけ相手にしたように書いてますが、実験はゲノム全体の遺伝子をターゲットにランダムかつ網羅的に行われており、mutant phenotypeの影響を調べた遺伝子数自体は１万よりずっと多くなります。実験結果をまとめたFitness Browserデータベースには、アノテーションがついた遺伝子も含めて実験結果を閲覧、比較できるようになっています。

8/6,8/7 誤字脱字修正 

2019 6/25　Preprint追記

2019 9/5　36生物に増加を確認

2020 初頭　37生物に増加したのを確認

2020 10/11 41生物 

2020 12/15 42生物

2022 1/12 46生物

2024 2/08 48生物

2024 04/19 49生物

　何千もの（Thousands of）バクテリアゲノムがシーケンシングされ、数百万のタンパク質の推定アミノ酸配列が明らかにされてきた。これらのタンパク質のわずかな割合しか実験的に研究されておらず、ほとんどのタンパク質の機能は、実験的に特徴付けられたタンパク質との類似性から予測されているだけである。しかし、バクテリアタンパク質の約3分の1には、この手法でアノテーションを付ける十分な相同性を持つ機能解析済みのタンパク質が存在しない（論文より ref.1）。さらに、これらの予測はしばしば不正確であり、これは相同なタンパク質が異なる基質特異性を有する可能性があることが１つの理由である（ref,2）。このsequence-to-function のギャップは微生物学にとってますます増大する課題である。なぜなら、新しいバクテリアゲノムは常に増加する速度でシーケンシングされているが、実験によるタンパク質の機能解析のスピードは比較的遅いままであるからである（ref.1）。

　未知のタンパク質の機能を研究するための1つのアプローチは、複数の条件下での対応する遺伝子の機能喪失突然変異の結果を評価することである（ref,3,4,5,6）。突然変異の表現型を比較ゲノミクスと組み合わせ、タンパク質の一部について証拠に基づきアノテーションを付与することができる（ref.3,4）。トランスポゾン突然変異誘発とその後のシークエンシング（TnSeq）は、何千何万もの異なる突然変異株が混ざったまま増殖する単一の実験からゲノムワイドの突然変異表現型を測定する（ref.7, 8）。 TnSeqと各変異株のrandom DNA barcoding（RB-TnSeq）とのカップリングは、多くの条件で表現型を測定することを容易にする（ref.9）。この論文では、RB-TnSeqを使用して、複数の実験条件下で32のバクテリアのそれぞれから、何千もの（Thousands of）遺伝子の突然変異表現型を系統的に探索し、配列間のギャップを解明する（論文より　図1a）。

• TnSeqについて（Tn-seq; high-throughput~より）

１、抽出したgenomic DNAにミニトランスポゾン（nptIIによるカナマイシン耐性あり）とトランスポザーゼを作用させ（認識サイトは"TA"）、ミニトランスポゾンカセット入りgenomic DNAを作る。

２、1で得た transposed DNAを使いバクテリアを形質転換する。

３、形質転換したバクテリアを様々な条件で継代する。

４、バクテリアから選抜前後のDNAを抽出する。

５、MmeI（リンク）（認識サイトから18-20 bp離れた部位を切る）で消化する。トランスポゾンカセットは両端に20 bp導入された部位のゲノム配列を付けた状態で切り出される（cの真ん中のイラスト）。

６、図のprimer1と2でPCRをかけ、プライマーとゲノムを含む領域をエンリッチする（cのイラストの赤い領域はマーカーであり、シーケンシングの必要はない）。

７、次世代にかけ、（アダプターを除き、）リファレンスゲノムにmappingする。20 bpしかないので、間違ってmappingされる可能性もあるが、そのような部位は実験から除外される。

８、指定条件でフィルタリングし（*1）、TMMなどで正規化後（ツール例）、遺伝子ごとに、培養前後でタグの数が増減したか調べる。

Tn Seqの概略。ref.8（Tn-seq; high-throughput~:  Tim van Opijnen, et al）より転載。

本論文では、TnSeqを何万回（バクテリア数  x 実験した条件）と行い、結果をまとめている。

• メモ1:　培養条件やgrowth mediaは、supplementaty tableファイル（S18）や、Fitness Browserの各バクテリア及び各条件のページに記載されている。本番の条件は、あらかじめ96-well（ 150µl volume）で培養を行い、TECAN infinite 200 PRO（リンク）などでOD600をモニターし、効果的かどうか調べた上で選択されている（*2）。形質転換方法（conjugationかelectroporation）はtable S20参照。大半の疑問に対する回答は論文のmethod　 sectionに書かれている。著者らの１つ前の論文も読んでおく（リンク）。

• メモ2

  例えばある遺伝子Aのトランスポゾンタグが、指定条件の培養開始前は100で、培養後は0になっているなら、遺伝子Aはその条件への適応に必須な遺伝子であることが強く示唆される。そもそも培養前から0なら、遺伝子Aは生育に必須なessential geneであることが強く示唆される。ある条件での培養後に有意に増加しているなら、遺伝子Aは指定した培養条件では増殖を抑える負の役割を持っていたことが示唆される（*4）。そのため、色々な部位にタグが入った各々異なるlocusのノックアウト株をごちゃまぜにして（poolして）しばらく培養すると、その遺伝子破壊によって分裂速度に変化が生じることで、培養後の遺伝子破壊株の数に差が生じる（いわゆるダーウィン進化が起きる）。言い換えると、locusにタグを持つセルの適応度に応じて遺伝子ごとに見つかるタグ数が増減する（あるlocusにタグが入った細胞はすごく減り、別のlocusにタグが入った細胞は増える）。このタグ数変化を、論文のFitness BrowserではFitness valuesという単位に正規化して変換し（0が変化なし）、遺伝子ごとにデータベース化している。 各条件、3万〜50万の変異株のプールから調べている。総数は実験条件とバクテリア数の掛け算になるので、数百条件を32バクテリアで行うと、とんでもない量の実験をこなさなければならない。Supplementaryデータをみると、データが細かすぎて目がおかしくなる。

実際にはバクテリアによって少し改変したプロトコルが用いられている。詳細はhttps://mbio.asm.org/content/6/3/e00306-15.full参照。

• Fitness Browser ヘルプ

http://fit.genomics.lbl.gov/cgi-bin/help.cgi

• 登録されているバクテリア:　40種以上のデータが利用可能（Organisms）。

シュードモナス属など８種出てきているが、これは、シュードモナスの多様性が大きく、種間、属間の遺伝子保存性が低いため、比較時の精度を上げる目的で使われている（論文より）。

使い方

Fitness Browser にアクセスする。

例えばαプロテオバクテリアのAzospirillum brasilense Sp245のデータを開いてみる。

nitrogen sourceが35と最も多い。nitrogen sourceの35をクリックしてみる。

使われた窒素源一覧（実験が成功したものだけ掲載されている）が表示される。

一番上の窒素源がAmmonium chlorideの実験をクリック。

画面が切り替わる。この条件で特異的な変化を示した遺伝子を調べるにはSpecificをクリック。

２遺伝子検出されている。上のAZOBR_RS25125をクリック。

fitnessが-2以上なのは、上２つのAmmonium chlorideだけなのが確認できる。

fitnessは、その条件でgrowthへの影響がなければゼロになる。0以下のfitness値はその条件へのfitnessが下がったことを意味し、0以上のfitness値は、その条件へのfitnessが上がっていることを意味する。ヘルプによれば、-1から1の範囲はばらつきで結果の解釈は不明瞭であるが、-2以下、または２以上になると強いfitness効果があると記載されている。値の信頼度は隣にt-like test（t値やZ scoreと同じスケールの値）で示される。データが数千あるため、tが-4以下、または4以上（すなわち|t| > 4）をsiginificantと記載している。

次に+ Geneタブをクリックする。200遺伝子以上検出されているが、fitnessが最大の遺伝子でもfitnessは1.8に留まっている（ただしt scoreは２を超えている）。

by conditionから調べても同じことができる。また、by Genesからは、遺伝子名で検索できる。

Fitness Browser - Exp Search

試しに、by Genesから、pyruvate dehydrogenase complexを調べてみる。生物種はさっきと同じAzospirillum brasilense Sp245とする。

E1ユニットをコードする遺伝子をクリック。

fitnessは "no data"になっている。これはTnseqでシーケンスされたタグがゼロまたはゼロ近くでデータがえられなかったか（=> essential gene）、実験条件を満たせなかったため、データから除外された可能性を表している（*1）。

近隣の遺伝子も表示されている。よく見ると、E1の上流には、画面に収まる範囲でno dataになっていないORFが３つある（しかし３つともphenotypeは"insignificant"）。その内、２つはhypothetical proteinとなっている。hypothetical proteinの１つAZOBR_RS22225の行の右端のinsigをクリックしてみる。

調べてみると、hypothetical proteinではあるが、炭素源、窒素源、ストレス（飢餓、抗生物質、温度変化などなど）、mobilityなど様々な条件でfitnessが増減していた。

近接のORFはどうなっているだろうか？（同じオペロンユニットで共転写されていたり、同じ転写因子の制御可にある可能性もある）。隣接遺伝子のfitness一覧を見るには、NearByのタブをクリックする。 

AZOBR_RS22225の下流側に同じ向きで隣接するAZOBR_RS22220は、fittingのパターンが同じとは言い難い。AZOBR_RS22225とAZOBR_RS22220はオペロンかもしれないが、少なくとも機能は異なる可能性が高まった。AZOBR_RS22225の上流側に同じ向きで隣接するAZOBR_RS2223については、Tnseqのデータが出ていない。AZOBR_RS2223をクリックするとNCBIなどから引っ張ってきているアノテーションの詳細が表示される。DNA topoisomerase IVのアノテーションがついているので、必須遺伝子のためtagのシーケンシングが出なかったのだろうか（ORFサイズが大きいこともessential geneの可能性が高めている => シーケンスデプスにもよるが、サイズが大きいほど、たまたまtag数がゼロか非常に少なかったという可能性は低くなる）...などなど、色々考えられる。

次に、距離を問わず同じfitness傾向を示すORFを探すため、Cofitのタブをクリックする。AZOBR_RS22225と同じ傾向を示すTop cofit リストが表示される。同じ傾向を示した遺伝子は、（著者らによれば）同じpathwayで機能する可能性が高い。

該当するバクテリアを扱っていない場合、BLAST解析して、Fitness Browserで調べられた遺伝子から、クエリとホモロジーが高い遺伝子をpick upできる（似た一次構造をもつタンパク質ならfitnessも似ている可能性が高い）。

BLAST

ORFのアミノ酸配列、または塩基配列をペーストしてsearchをクリック。

しばらくするとblast結果が出てくる。e-valueの閾値を超えたリストが並ぶ。

 十分な相同性があるなら、同様のmutant phenotypeを示す可能性も高くなる。

他にもいくつかの機能があり、tophitを抽出してきたり、選択したデータから図を描いたりできます。

感想

ものすごい実験の数とデータ量に圧倒されますが、その分、使いこなせば、武器になるデータベースです。（なんせ、個々のバクテリアの研究者が実験を始める前に論文の著者らがknock out株を作ってgrowthが変化するかどうかをいろんな条件で手当たり次第に調べ、さらに、そのデータを比較可能な形で公開してくれているわけです。とんでもない話です）。

このデータベースの使い方ですが、いきなりhypothetical proteinを探し始めるより、まずはこれまで、または現在進行形で機能解析している遺伝子や、第三者によって機能がよく調べられきた遺伝子を使って、どんなfitnessを示すか調べてみて下さい。遊んでいるうちに使い方を学べることもありますが、機能解析が幅白い文脈で調べられている遺伝子を使う方が、この手法の可能性と限界が見えてくるんじゃないかと思います(*5)。

引用

Mutant phenotypes for thousands of bacterial genes of unknown function

Price MN, Wetmore KM, Waters RJ, Callaghan M, Ray J, Liu H, Kuehl JV, Melnyk RA, Lamson JS, Suh Y, Carlson HK, Esquivel Z, Sadeeshkumar H, Chakraborty R, Zane GM, Rubin BE, Wall JD, Visel A, Bristow J, Blow MJ, Arkin AP, Deutschbauer AM

Nature. 2018 May;557(7706):503-509. 

https://www.natureasia.com/ja-jp/nature/557/7706/s41586-018-0124-0/機能が未知の数千の細菌遺伝子の変異表現型

*1

論文の"Identifying essential or nearly essential genes"のパラグラフを参照。タグの分布をPoisson distributionと仮定し、遺伝子のサイズから偶然tagがゼロになる可能性が棄却された遺伝子のみ、Fitness Browserにデータが示されている。

*2

補足表には、本実験以上を断念した条件などもまとめられている（S4）。例えば、カナマイシンにnativeに耐性があった、など。

*3

Tn-Seqを説明した動画:

Tn-seq in Rhodobacter sphaeroides - mSystems®（American Society for Microbiology）

          2:20ごろに実際の挿入状況。 

*4 

変動するのは生物学的要因以外も考えられる。バイアスは先行研究で議論されている（pubmed）

追記

*5

in vitroの解析で機能がわかっていると思っていたタンパク質に別の機能も備わっていることがこのデータから見えてくるかもしれません。

2019 6/25 追記

Large-scale chemical-genetics of the human gut bacterium Bacteroides thetaiotaomicron

Hualan Liu, Morgan N. Price, Hans K. Carlson, Yan Chen, Jayashree Ray, Anthony L. Shiver, Christopher J. Petzold, Kerwyn Casey Huang, Adam P. Arkin, Adam M. Deutschbauer
doi: https://doi.org/10.1101/573055

のデータが追記され、現在35種のデータを利用できるようになっています。

関連

GapMindやKEGGなどへのリンクも追加されています。

sam/bamがmalformedではないか調べるPicardのValidateSamFile

sam/bamをいじっていると、ヘッダーが無かったり重複したり、ダウンロードが不完全だったり、様々な理由でおかしくなってしまうことがある。PicardのValidateSamFileはsam/bamにエラーがないか分析するコマンド。実行するとエラーが見つかったところを教えてくれる。

GATKより。

https://software.broadinstitute.org/gatk/documentation/article.php?id=7571

Picard-toolsのインストール

 brewで導入できる。

brew install picard-tools

> picard ValidateSamFile

$ picard ValidateSamFile

ERROR: Option 'INPUT' is required.

USAGE: ValidateSamFile [options]

Documentation: http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

Validates a SAM or BAM file.

This tool reports on the validity of a SAM or BAM file relative to the SAM format specification.  This is useful for

troubleshooting errors encountered with other tools that may be caused by improper formatting, faulty alignments,

incorrect flag values, etc. 

By default, the tool runs in VERBOSE mode and will exit after finding 100 errors and output them to the console

(stdout). Therefore, it is often more practical to run this tool initially using the MODE=SUMMARY option.  This mode

outputs a summary table listing the numbers of all 'errors' and 'warnings'.

When fixing errors in your file, it is often useful to prioritize the severe validation errors and ignore the

errors/warnings of lesser concern.  This can be done using the IGNORE and/or IGNORE_WARNINGS arguments.  For helpful

suggestions on error prioritization, please follow this link to obtain additional documentation on ValidateSamFile

(https://www.broadinstitute.org/gatk/guide/article?id=7571).

After identifying and fixing your 'warnings/errors', we recommend that you rerun this tool to validate your SAM/BAM file

prior to proceeding with your downstream analysis.  This will verify that all problems in your file have been addressed.

Usage example:

java -jar picard.jar ValidateSamFile \

I=input.bam \

MODE=SUMMARY

To obtain a complete list with descriptions of both 'ERROR' and 'WARNING' messages, please see our additional 

documentation (https://www.broadinstitute.org/gatk/guide/article?id=7571) for this tool.

Version: 2.18.9-SNAPSHOT

Options:

--help

-h                            Displays options specific to this tool.

--stdhelp

-H                            Displays options specific to this tool AND options common to all Picard command line

                              tools.

--version                     Displays program version.

INPUT=File

I=File                        Input SAM/BAM file  Required. 

OUTPUT=File

O=File                        Output file or standard out if missing  Default value: null. 

MODE=Mode

M=Mode                        Mode of output  Default value: VERBOSE. This option can be set to 'null' to clear the

                              default value. Possible values: {VERBOSE, SUMMARY} 

IGNORE=Type                   List of validation error types to ignore.  Default value: null. Possible values:

                              {INVALID_QUALITY_FORMAT, INVALID_FLAG_PROPER_PAIR, INVALID_FLAG_MATE_UNMAPPED,

                              MISMATCH_FLAG_MATE_UNMAPPED, INVALID_FLAG_MATE_NEG_STRAND, MISMATCH_FLAG_MATE_NEG_STRAND,

                              INVALID_FLAG_FIRST_OF_PAIR, INVALID_FLAG_SECOND_OF_PAIR,

                              PAIRED_READ_NOT_MARKED_AS_FIRST_OR_SECOND, INVALID_FLAG_NOT_PRIM_ALIGNMENT,

                              INVALID_FLAG_SUPPLEMENTARY_ALIGNMENT, INVALID_FLAG_READ_UNMAPPED, INVALID_INSERT_SIZE,

                              INVALID_MAPPING_QUALITY, INVALID_CIGAR, ADJACENT_INDEL_IN_CIGAR, INVALID_MATE_REF_INDEX,

                              MISMATCH_MATE_REF_INDEX, INVALID_REFERENCE_INDEX, INVALID_ALIGNMENT_START,

                              MISMATCH_MATE_ALIGNMENT_START, MATE_FIELD_MISMATCH, INVALID_TAG_NM, MISSING_TAG_NM,

                              MISSING_HEADER, MISSING_SEQUENCE_DICTIONARY, MISSING_READ_GROUP, RECORD_OUT_OF_ORDER,

                              READ_GROUP_NOT_FOUND, RECORD_MISSING_READ_GROUP, INVALID_INDEXING_BIN,

                              MISSING_VERSION_NUMBER, INVALID_VERSION_NUMBER, TRUNCATED_FILE,

                              MISMATCH_READ_LENGTH_AND_QUALS_LENGTH, EMPTY_READ, CIGAR_MAPS_OFF_REFERENCE,

                              MISMATCH_READ_LENGTH_AND_E2_LENGTH, MISMATCH_READ_LENGTH_AND_U2_LENGTH,

                              E2_BASE_EQUALS_PRIMARY_BASE, BAM_FILE_MISSING_TERMINATOR_BLOCK, UNRECOGNIZED_HEADER_TYPE,

                              POORLY_FORMATTED_HEADER_TAG, HEADER_TAG_MULTIPLY_DEFINED,

                              HEADER_RECORD_MISSING_REQUIRED_TAG, HEADER_TAG_NON_CONFORMING_VALUE, INVALID_DATE_STRING,

                              TAG_VALUE_TOO_LARGE, INVALID_INDEX_FILE_POINTER, INVALID_PREDICTED_MEDIAN_INSERT_SIZE,

                              DUPLICATE_READ_GROUP_ID, MISSING_PLATFORM_VALUE, INVALID_PLATFORM_VALUE,

                              DUPLICATE_PROGRAM_GROUP_ID, MATE_NOT_FOUND, MATES_ARE_SAME_END,

                              MISMATCH_MATE_CIGAR_STRING, MATE_CIGAR_STRING_INVALID_PRESENCE,

                              INVALID_UNPAIRED_MATE_REFERENCE, INVALID_UNALIGNED_MATE_START, MISMATCH_CIGAR_SEQ_LENGTH,

                              MISMATCH_SEQ_QUAL_LENGTH, MISMATCH_FILE_SEQ_DICT, QUALITY_NOT_STORED, DUPLICATE_SAM_TAG,

                              CG_TAG_FOUND_IN_ATTRIBUTES} This option may be specified 0 or more times. 

MAX_OUTPUT=Integer

MO=Integer                    The maximum number of lines output in verbose mode  Default value: 100. This option can be

                              set to 'null' to clear the default value. 

IGNORE_WARNINGS=Boolean       If true, only report errors and ignore warnings.  Default value: false. This option can be

                              set to 'null' to clear the default value. Possible values: {true, false} 

VALIDATE_INDEX=Boolean        DEPRECATED.  Use INDEX_VALIDATION_STRINGENCY instead.  If true and input is a BAM file

                              with an index file, also validates the index.  Until this parameter is retired VALIDATE

                              INDEX and INDEX_VALIDATION_STRINGENCY must agree on whether to validate the index. 

                              Default value: true. This option can be set to 'null' to clear the default value. Possible

                              values: {true, false} 

INDEX_VALIDATION_STRINGENCY=IndexValidationStringency

                              If set to anything other than IndexValidationStringency.NONE and input is a BAM file with

                              an index file, also validates the index at the specified stringency. Until VALIDATE_INDEX

                              is retired, VALIDATE INDEX and INDEX_VALIDATION_STRINGENCY must agree on whether to

                              validate the index.  Default value: EXHAUSTIVE. This option can be set to 'null' to clear

                              the default value. Possible values: {EXHAUSTIVE, LESS_EXHAUSTIVE, NONE} 

IS_BISULFITE_SEQUENCED=Boolean

BISULFITE=Boolean             Whether the SAM or BAM file consists of bisulfite sequenced reads. If so, C->T is not

                              counted as an error in computing the value of the NM tag.  Default value: false. This

                              option can be set to 'null' to clear the default value. Possible values: {true, false} 

MAX_OPEN_TEMP_FILES=Integer   Relevant for a coordinate-sorted file containing read pairs only. Maximum number of file

                              handles to keep open when spilling mate info to disk. Set this number a little lower than

                              the per-process maximum number of file that may be open. This number can be found by

                              executing the 'ulimit -n' command on a Unix system.  Default value: 8000. This option can

                              be set to 'null' to clear the default value. 

SKIP_MATE_VALIDATION=Boolean

SMV=Boolean                   If true, this tool will not attempt to validate mate information. In general cases, this

                              option should not be used.  However, in cases where samples have very high duplication or

                              chimerism rates (> 10%), the mate validation process often requires extremely large

                              amounts of memory to run, so this flag allows you to forego that check.  Default value:

                              false. This option can be set to 'null' to clear the default value. Possible values:

                              {true, false} 

ラン 

bamを分析する。MODE=SUMMARYを外すと、異常なリードを全てプリントする。

picard ValidateSamFile I=input.bam MODE=SUMMARY 

$ picard ValidateSamFile I=~/Documents/input.sam MODE=SUMMARY

#一部略

MAX_OPEN_TEMP_FILES=8000 SKIP_MATE_VALIDATION=false VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json USE_JDK_DEFLATER=false USE_JDK_INFLATER=false

[Sat Aug 04 17:26:55 JST 2018] Executing as kamisakakazuma@kamisakBookpuro on Mac OS X 10.13.6 x86_64; OpenJDK 64-Bit Server VM 1.8.0_92-b15; Deflater: Intel; Inflater: Intel; Provider GCS is not available; Picard version: 2.18.9-SNAPSHOT

No errors found

sam/bamに異常が見つからなければ、最後に"No errors found"が表示される。問題がある場合、WARNINGSかERRORSで内容が示される。

bamのリードグループを改変してプログラムの挙動を確かめる。

samtools view -H input.bam > modify_header.sam
vi modify_header.sam

オリジナル

@RG     ID:X    LB:Y    SM:sample1      PL:ILLUMINA

改変後

@RG     ID:X    LB:Y    SM:sample1      SM:sample1      PL:ILLUMINA

レアケース過ぎて例としては良くないが、手っ取り早くサンプル名を2つにしてみた。

修正ヘッダを元のbamヘッダー情報と置換する。

samtools reheader modify_header sinput.bam > reheaer.bam

ValidateSamFileを実行。

picard ValidateSamFile I=reheaer.bam MODE=SUMMARY 

ERROR:HEADER_TAG_MULTIPLY_DEFINED内容も含めてエラーを検出できている。

個別のケースの修正方法については、GAATKチュートリアルで書かれている例を参考にして下さい。

GATK | Doc #7571 | Errors in SAM/BAM files can be diagnosed with ValidateSamFile

引用

Picard Tools - By Broad Institute

アンプリコンシーケンシングのdenoiseを行う DUDE-Seq

　次世代シーケンシング（NGS）技術と呼ばれる新しい世代の高スループット、低コストのシークエンシング技術が、大規模な比較研究や進化研究などの生物医学研究を再構成している[ref.2-4]。自動サンガーシーケンシングと比較すると、NGSプラットフォームは、大幅に短いリードが大量に生成され、これは新しいさまざまな計算上の課題を引き起こす[ref.5]。

　全ゲノムシーケンシング（WGS）、クロマチン免疫沈降（ChIP）シーケンシング、およびターゲットシーケンシングを含むNGS [ref.6,7]を用いるいくつかのDNAシーケンシング法がある。 WGSは、生物のゲノムを解析してすべてのバリアントを捕獲し、潜在的な原因変異を同定するために使用される。De novoゲノムアセンブリにも使用されている。ChIPシーケンシングは、転写因子および他のタンパク質のゲノムワイドDNA結合部位を同定する。この論文の焦点であるターゲットシーケンシング（exome sequencing and amplicon sequencing）は、研究者が特定の表現型に関与する可能性のある関心領域を調べることに焦点を当てることができる費用対効果の高い方法である。以前の研究[ref.8,9]によれば、ターゲットシーケンシングはしばしば疾患関連遺伝子のエキソンの完全なカバレッジをもたらし、代替の方法は約90〜95％のカバレッジをもたらす。したがって、臨床の現場では、研究者は診断評価のためターゲットシーケンシングに頼っている傾向がある。

　分子レベルで蛍光標識に基づいて配列を検出するため、NGS技術は、通常、エマルジョンPCRまたはsolid-phase amplification[ref.1]によって増幅される鋳型を必要とするイメージングシステムに依存する。これらの増幅およびイメージングプロセスはエラーの入ったリード生成する。それが起源となって、誤ったホモポリマー長、挿入/欠失/置換の導入、およびPCRキメラ[ref.6]が生じ得る。イルミナを含む多くのプラットフォームでは置換エラーが支配的であり、挿入および欠失（indels）によるホモポリマーエラーは、454パイロシーケンシングおよびイオントレントでも豊富である。

　エラーを含むリードは、下流の分析（例えば、バリアントコールおよびゲノムアセンブリ）を複雑にし、解析パイプライン全体の品質を低下させるため、適切に処理する必要がある。Soft-clipping[ref.7]はリードの各塩基のクオリティに応じてリードの3 '末端をトリムする（クオリティコントロールの）最も単純なアプローチだが、情報が失われる[ref.10]。より洗練された方法は、シーケンスデータのエラーの検出と訂正に焦点を当てている[ref.11-20]。イルミナシーケンシングプラットフォームが広く使用されていることを考えると、ほとんどのエラー訂正アルゴリズムは置換エラーをターゲットにしている[ref.10]。

　最近のレビュー[ref.10,21]に要約されているように、NGSデータの現在のエラー訂正方法は以下のように分類できる: k-mer（すなわち、長さkのオリゴヌクレオチド）frequency/スペクトラムベース、multiple sequence alignment（MSA）ベース、および統計的な誤差モデルに基づく方法。 k-merベースの方法の背後にあるアイデアは、入力スペクトラムから信頼できるk-merリストを作成し、このスペクトラムによって表されるコンセンサスに基づいて信頼できないk-mersを修正することである[ref.13,20,22-25]。 k-merの長さに加えて、カバレッジ（k-mer発生）情報は、信頼できるk-mersを決定するために重要である。エラーが稀でランダムであり、カバレッジが均一であるという仮定の下で、十分に大きなkについては、ほとんどのエラーがk-mersをゲノム中に存在しないk-mer配列に変えることを期待することは合理的である。したがって、NGSによって得られた高カバレッジゲノム配列について、疑わしいk-merを同定し、コンセンサスに基づいてそれらを修正することができる。 MSAベースの方法[ref.12,16,26]は、関連する配列をそれらの類似性に従ってアライメントし、様々な技法、通常はアラインメントカラムのコンセンサスに基づいてアライメントしたリードを修正することで機能する。このアライメントベースのスキームは、indelエラーを修正するために本質的に適している。初期の方法は計算上の問題を抱えていたが、最近のアプローチでは高度なindexing手法を使用してアライメントを高速化している。統計的エラーモデルに基づく方法[ref.27-29]では、エラー発生を含んだシーケンシングプロセスを捕捉する統計的モデルが開発されている。これに関して、経験的な混同モデルは、例えば、アライメント結果、Phredクオリティスコア（自動DNAシーケンシングによって生成された核酸塩基のクオリティ尺度）[ref.30]から得られた情報、または他のパラメータを利用して、データセットから作成されることが多い。

　上記の方法は、多くの場合、さまざまなプラットフォームで良好なパフォーマンスを発揮するが、いくつかの制限もある。第1に、k-merベースのスキームは、トランスクリプトミックス、メタゲノミクス、ヘテロジニアスな細胞サンプル、または事前増幅されたライブラリ[ref.21]のように、カバレッジが変化すると予想される場合、不適格となる傾向がある。。第2に、不均一なカバレッジに関連する上記の問題を抱えていないMSAベースの方法は、アライメントした列に対するヒューリスティックかつ洗練されたコンセンサス決定規則の適用を必要とし、このような規則はしばしば特定のアプリケーションまたはシーケンスプラットフォームではセンシティブになるかもしれない。第3に、統計的なエラーモデルベースの方法は、基礎となるDNA配列に対する追加の確率論的モデリング仮定のために、計算上高価なスキーム（例えば、期待値最大化(EM)）を使用する。さらに、このようなモデリングの前提条件に対して妥当性と精度にはほとんど注意が払われない。さらに、ほぼ最適または健全なエラー訂正性能が達成されたかどうかの理論的分析が行われているかに対しては、言うまでもない。最後に、3つの方法を適用する多くの既存のスキームは、入力シーケンスをマージすることによって作成された代表のコンセンサスノイズ除去シーケンスのみを返すことが多い。したがって、消滅後の配列の数はしばしば保存されない。一部のアプリケーションでは、これによってダウンストリーム分析に矛盾が生じることがある。これらの限界に対処するために、多くの既存のツールがパフォーマンスを改善するために3つの方法を補完的に組み合わせている[ref.10、21]。

　この論文では、代わりに、正確なDNA配列のノイズ除去のためにDiscrete Universal DEnoiser（DUDE）[ref.31 リンク]というアルゴリズムを適用した。 DUDEは、各ソースシンボルを独立して統計的に同一に破壊するノイズメカニズムによって破損した有限値成分（ソースシンボル）を有する再構成シーケンスの一般的なセッティングのために開発された。 DNA denoising の文献では、そのようなノイズモデルは、統計的エラーモデルベースの方法で一般的に使用されるconfusion matrixと同等である。オリジナルの論文[ref.31]に示されているように、DUDEは以下の設定に対して厳しい性能保証を示している。既知のノイズメカニズムを前提とした場合にのみ、元のクリーンソースデータに対して確率的モデリング仮定が行われていない場合でも、データが増加するあらゆるソースデータに対して最適なノイズ除去性能を普遍的に達成することが示される。上記のDUDEの設定は、DNA配列のノイズ除去の設定に当てはまることに注意する。すなわち、クリーンなDNA配列の正確な確率モデルを確立するのは難しいが、リファレンス配列に基づいてシーケンシングデバイスのノイズモデル（confusion matrix）を仮定するのは簡単である。

（一段落省略）

上記のDUDEアルゴリズムのユニークな性質を利用して、著者らの実験では、他の最先端の手法よりも優れていることを示している。具体的には、ターゲットアンプリコンシーケンシング（例えば、ガン遺伝子、16S rRNA、植物および動物のシーケンシング[ref.32]）の適用可能な領域の中で、異なるライブラリー調製法およびDNAポリメラーゼで得られた16S rRNAベンチマークデータセットを使用して、アルゴリズムを使用して、ターゲットアンプリコンシーケンシングデータセットは、しばしばWGSやChIPデータセットよりも深いシーケンシングカバレッジを持っており、従来のk-merベースの手法では頻繁に増幅バイアス問題が発生する[ref.33]。対照的に、DUDE-Seqでは、シーケンスカバレッジが深くなるにつれて、コンテキストカウントベクトルがより確率の高いコンテキストを蓄積することがあり、ノイズ除去の堅牢性が向上する。 DUDEの2つのバージョンを、別々に、置換エラーとホモポリマーエラーに適用する（以下省略）。

公式サイト

http://data.snu.ac.kr/pub/dude-seq/

Github

 DUDE-Seqに関するツイート。

使い方

ここではweb版を紹介する。

fasta、fastq、またdatの圧縮ファイル（zipかgzip）を指定してsubmitする。指定できるパラメータはアルゴリズム1 (substitution-error correction) のkサイズとアルゴリズム２ (homopolymer-error correction)のkサイズ。

ジョブが終わると、denoiseされた配列がダウンロードできるようになる。

追記

ローカル版のインストール

ubuntu18.04のpython2.7.14でテストした。

依存

• python2
• libboost-dev
• libgsl0-dev
• liblapack-dev
• zlib1g-dev

pythonライブラリ

• biom-format（リンク）
• cogent
• numpy

sudo apt install libboost-dev libgsl0-dev liblapack-dev zlib1g-dev
pip install biom-format==1.3.0 #最新の2.1.6だとbiom.tableをimportできなかった
pip install cogent
pip install numpy
pip install qcli

本体　GIthub

git clone https://github.com/datasnu/dude-seq.git
cd dude-seq/

> python Scripts/process_sff.py -h

$ python Scripts/process_sff.py -h

Usage: process_sff.py [options] {-i/--input_dir INPUT_DIR}

[] indicates optional input (order unimportant)

{} indicates required input (order unimportant)

This script converts a directory of sff files into FASTA, QUAL and flowgram files.

Example usage: 

Print help message and exit

 process_sff.py -h

Simple example: Convert all the sffs in directory "sffs/" to fasta and qual.

 process_sff.py -i sffs/

Convert a single sff to fasta and qual.

 process_sff.py -i sffs/test.sff

Flowgram example: Convert all the sffs in directory "sffs/" to fasta and qual, along with a flowgram file.

 process_sff.py -i sffs/ -f

Convert a single sff to fasta and qual, along with a flowgram file.

 process_sff.py -i sffs/test.sff -f

Output example: Convert all the sffs in directory "sffs/" to fasta and qual, along with a flowgram file and write them to another directory.

 process_sff.py -i sffs/ -f -o output_dir

Options:

  --version             show program's version number and exit

  -h, --help            show this help message and exit

  -v, --verbose         Print information during execution -- useful for

                        debugging [default: False]

  --no_trim             do not trim sequence/qual (requires --use_sfftools

                        option) [default: False]

  -f, --make_flowgram   generate a flowgram file. [default: False]

  -t, --convert_to_FLX  convert Titanium reads to FLX length. [default: False]

  --use_sfftools        use the external programs sfffile and sffinfo for

                        processing, instead of the equivalent python

                        implementation

  -o OUTPUT_DIR, --output_dir=OUTPUT_DIR

                        Input directory of sff files [default: same as input

                        dir]

  REQUIRED options:

    The following options must be provided under all circumstances.

    -i INPUT_DIR, --input_dir=INPUT_DIR

                        Input directory of sff files or a single sff filepath

                        [REQUIRED]

追記

ローカル版のインストール

ubuntu18.04のpython2.7.14でテストした。

依存

• python2
• libboost-dev
• libgsl0-dev
• liblapack-dev
• zlib1g-dev

pythonライブラリ

• biom-format（リンク）
• cogent
• numpy

sudo apt install libboost-dev libgsl0-dev liblapack-dev zlib1g-dev
pip install biom-format==1.3.0 #最新の2.1.6だとbiom.tableをimportできなかった
pip install cogent
pip install numpy
pip install qcli

本体　GIthub

git clone https://github.com/datasnu/dude-seq.git
cd dude-seq/

> python Scripts/process_sff.py -h

$ python Scripts/process_sff.py -h

Usage: process_sff.py [options] {-i/--input_dir INPUT_DIR}

[] indicates optional input (order unimportant)

{} indicates required input (order unimportant)

This script converts a directory of sff files into FASTA, QUAL and flowgram files.

Example usage: 

Print help message and exit

 process_sff.py -h

Simple example: Convert all the sffs in directory "sffs/" to fasta and qual.

 process_sff.py -i sffs/

Convert a single sff to fasta and qual.

 process_sff.py -i sffs/test.sff

Flowgram example: Convert all the sffs in directory "sffs/" to fasta and qual, along with a flowgram file.

 process_sff.py -i sffs/ -f

Convert a single sff to fasta and qual, along with a flowgram file.

 process_sff.py -i sffs/test.sff -f

Output example: Convert all the sffs in directory "sffs/" to fasta and qual, along with a flowgram file and write them to another directory.

 process_sff.py -i sffs/ -f -o output_dir

Options:

  --version             show program's version number and exit

  -h, --help            show this help message and exit

  -v, --verbose         Print information during execution -- useful for

                        debugging [default: False]

  --no_trim             do not trim sequence/qual (requires --use_sfftools

                        option) [default: False]

  -f, --make_flowgram   generate a flowgram file. [default: False]

  -t, --convert_to_FLX  convert Titanium reads to FLX length. [default: False]

  --use_sfftools        use the external programs sfffile and sffinfo for

                        processing, instead of the equivalent python

                        implementation

  -o OUTPUT_DIR, --output_dir=OUTPUT_DIR

                        Input directory of sff files [default: same as input

                        dir]

  REQUIRED options:

    The following options must be provided under all circumstances.

    -i INPUT_DIR, --input_dir=INPUT_DIR

                        Input directory of sff files or a single sff filepath

                        [REQUIRED]

引用

DUDE-Seq: Fast, flexible, and robust denoising for targeted amplicon sequencing
Byunghan Lee, Taesup Moon, Sungroh Yoon, Tsachy Weissman

PLoS One. 2017; 12(7): e0181463. Published online 2017 Jul 27.

複数のSV検出結果を統合し、精度の高いSVコールを行う MetaSV

8/19 pindel、lumpy、cnvkit、breakdanerコマンドミス修正

　SVは、ゲノムの多様性およびゲノムの障害に寄与することに関与している（Stankiewicz and Lupski、2010）。したがって、SVの検出には相当量の作業が行われている。一般に、SVを検出するためのツールは、リードのアライメントからの以下の信号のうちの1つまたは複数を使用する：split-read [分割アライメント。例: Pindel（Ye et、2009）]、read-pair[異常なペアエンドアライメント。例: BreakDancer（Chen et al、2009）]、カバレッジデプス[異常なカバレージ。例: CNVnator（Abyzov et al、2011）]、ジャンクションマッピング[既知のSVブレークポイントへのアライメント。例: BreakSeq2（Abyzov et al、2015; Lam et al、2010）]またはブレークポイント周辺のアセンブリ[MindTheGap（Rizk et al、2014）]。しかし、SVタイプとサイズのニッチがあり、すべてのタイプのSVを包括的に検出するシグナルはない。これは、SV検出を改善するための複数の方法を統合するツールの開発を必要とする。

　以前の研究（Lam et al、2012; Mills et al、2011）は、複数のツールと方法によって作成されたバリアントコールが一般により正確であることを示している。このため、複数の方法を採用するためのツールが開発されている。 DELLY（Rausch et al、2012）、LUMPY（Layer etal、2014）、SVMerge（Wong etal、2010）などの複数のツールの結果をマージするものなどがある。ただし、LUMPYとDELLYは挿入を検出することができず、SVMergeはコールのマージの際に個々のツールのSV解像度を無視する。したがって、この研究は、異なるタイプとサイズのSVを高い精度と解像度で検出するための既存のSVマージングツールの限界に取り組もうとしている。

　MetaSVは、複数のSV検出方法とツールを使用して、信頼性が高く正確なSVコールセットを検出する。 MetaSVの新しさは、以下の重要なアイデアの組み合わせにある。複数の独立した方法の組み合わせで報告されたコールは一般に品質が良く、動的プログラミングによるローカルアセンブリを使用してSVブレークポイントを絞り込むことができる。

　MetaSVは、以下のステップ（論文　図1）で進行する。

Intra-tool merging：同じツールで生成された潜在的なduplicate callsがここでマージされる。十分なオーバーラップがある場合、2つのコールはマージされる。
Inter-tool merging：複数のツールで生成されたコールがマージされる。複数のツールに共通するコールのブレークポイントを決定する際には、精度が高いことが知られているメソッドが優先される（e.g., スプリット・リード情報の方がリード・ペア情報より優先）。このメソッド認識マージは、報告されたSVのブレークポイントが可能な限り正確であることを確実にするために重要である。
Local assembly：追加の証拠を収集してSVの配列を決定するため、ツールによって報告されたSV領域に対してローカルアセンブリが実行される。アセンブリのため独自のペア情報を使用する機能があるため、SPAdes（Bankevich etal、2012）アセンブラが使用される。
Breakpoint resolution：アセンブリされたSV配列は、ダイナミックプログラミングを使用してブレークポイントを検出または調整するためリファレンスとアライメントされる（Abyzov and Gerstein、2011）。
Genotyping：SVブレークポイント周辺のカバレッジがSVのzygosityを決定するために使用される。最終的な出力は、各SVのgenotypeを持つVCFファイルとして生成される。
Annotation：MetaSVは入力と出力をVCFで標準化している。各SVには、個々のツールによって行われた対応するコーラーを示し、その信頼レベルを分類するために注釈が付けられる。複数のツールで検出されたSVは、高い信頼性があるとみなされる。

MetaSV HP

http://bioinform.github.io/metasv/

MetaSV: An Accurate and Integrative Structural-Variant Caller for Next Generation Sequencing

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

 MetaSVに関するツイート。

インストール

mac os 10.12のAnaconda環境でテストしたが、アセンブリをonにするとエラーを起こした（spades 3.6.2と3.12でテストした）。

依存

• pybedtools: Note that bedtools (version 2.21 or greater) has to be installed separately in order for pybedtools to work. Please use pybedtools version 0.6.9.
• pysam-0.7.7: MetaSV has only been tested with version 0.7.7 of pysam.
• pyvcf: Please use pyvcf version 0.6.7.

In addition, paths to the following tools must be provided as MetaSV arguments:

• SPAdes: Used for assembly around the SV breakpoints
• AGE: Used for determining the SV breakpoints using assembled contigs. Please compile AGE without OpenMP support using make OMP=no.

spades（bioconda）、age（bioconda）もcondaで導入可能（"spades.py"と "age_align"でヘルプが出るなら、whichでパス確認）。

spadesは"conda search spades"でインストール可能なバージョンを調べ、3.6.2を導入した(*5)。

conda search spades
conda install -c bioconda spades==3.6.2

上記に加え、metaSVに入力するSVをコールするツールも必要になる。 ヘルプには以下のSV検出ツール（SV caller）のオプションが記載されている。

• Pindel
• breakdancer (breakdancer-max)
• breakseq (2も含む?)
• cnvnator
• gatk
• manta
• lumpy
• cnvkit
• wham

SV検出に使うcallerだけ導入すればよい。インストール方法については、各コマンド実行部分に簡単な説明を残した。

本体 Github

condaで導入できる（リンク）。

conda install -c bioconda metasv 
#permittionエラーが出たので、てテスト時はsudoで0.4.0を入れた (" metasv==0.4.0")

>  run_metasv.py

$ run_metasv.py -h

usage: run_metasv.py [-h] --sample Sample

                     [--pindel_vcf pindel_vcf [pindel_vcf ...]]

                     [--pindel_native File list [File list ...]]

                     [--breakdancer_vcf breakdancer_vcf [breakdancer_vcf ...]]

                     [--breakdancer_native File list [File list ...]]

                     [--breakseq_vcf breakseq_vcf [breakseq_vcf ...]]

                     [--breakseq_native breakseq_native [breakseq_native ...]]

                     [--cnvnator_vcf cnvnator_vcf [cnvnator_vcf ...]]

                     [--cnvnator_native File list [File list ...]]

                     [--gatk_vcf file [file ...]]

                     [--manta_vcf MANTA_VCF [MANTA_VCF ...]]

                     [--lumpy_vcf LUMPY_VCF [LUMPY_VCF ...]]

                     [--cnvkit_vcf CNVKIT_VCF [CNVKIT_VCF ...]]

                     [--wham_vcf WHAM_VCF [WHAM_VCF ...]]

                     [--mean_read_length MEAN_READ_LENGTH] --reference

                     reference [--chromosomes CHROMOSOMES [CHROMOSOMES ...]]

                     [--gaps gaps] [--filter_gaps] [--keep_standard_contigs]

                     [--bams BAMS [BAMS ...]] [--isize_mean ISIZE_MEAN]

                     [--isize_sd ISIZE_SD] [--wiggle WIGGLE]

                     [--inswiggle INSWIGGLE] [--minsvlen MINSVLEN]

                     [--maxsvlen MAXSVLEN] [--overlap_ratio OVERLAP_RATIO]

                     [--min_avg_base_qual MIN_AVG_BASE_QUAL]

                     [--min_mapq MIN_MAPQ] [--min_soft_clip MIN_SOFT_CLIP]

                     [--max_nm MAX_NM] [--min_matches MIN_MATCHES]

                     [--min_support_ins MIN_SUPPORT_INS]

                     [--min_support_frac_ins MIN_SUPPORT_FRAC_INS]

                     [--max_ins_intervals MAX_INS_INTERVALS]

                     [--mean_read_coverage MEAN_READ_COVERAGE]

                     [--min_ins_cov_frac MIN_INS_COV_FRAC]

                     [--max_ins_cov_frac MAX_INS_COV_FRAC]

                     [--sc_other_scale SC_OTHER_SCALE] [--spades SPADES]

                     [--spades_options SPADES_OPTIONS]

                     [--spades_timeout SPADES_TIMEOUT] [--disable_assembly]

                     [--svs_to_assemble {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]]

                     [--svs_to_softclip {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]]

                     [--extraction_max_read_pairs EXTRACTION_MAX_READ_PAIRS]

                     [--spades_max_interval_size SPADES_MAX_INTERVAL_SIZE]

                     [--assembly_max_tools ASSEMBLY_MAX_TOOLS]

                     [--assembly_pad ASSEMBLY_PAD] [--stop_spades_on_fail]

                     [--age AGE] [--age_timeout AGE_TIMEOUT]

                     [--min_inv_subalign_len MIN_INV_SUBALIGN_LEN]

                     [--min_del_subalign_len MIN_DEL_SUBALIGN_LEN]

                     [--age_window AGE_WINDOW] [--boost_sc]

                     [--gt_window GT_WINDOW] [--gt_normal_frac GT_NORMAL_FRAC]

                     [--svs_to_report {INV,CTX,INS,DEL,ITX,DUP} [{INV,CTX,INS,DEL,ITX,DUP} ...]]

                     [--enable_per_tool_output] [--workdir WORKDIR]

                     [--num_threads NUM_THREADS] --outdir OUTDIR [--version]

Merge SVs from multiple tools for accurate SV calling

optional arguments:

  -h, --help            show this help message and exit

Input data options:

  --sample Sample       Sample name (default: None)

  --pindel_vcf pindel_vcf [pindel_vcf ...]

                        VCF file or dir for Pindel VCFs (default: )

  --pindel_native File list [File list ...]

                        Pindel native files (default: )

  --breakdancer_vcf breakdancer_vcf [breakdancer_vcf ...]

                        VCF file or dir for BreakDancer VCFs (default: )

  --breakdancer_native File list [File list ...]

                        BreakDancer native files (default: )

  --breakseq_vcf breakseq_vcf [breakseq_vcf ...]

                        VCF file or dir for BreakSeq VCFs (default: )

  --breakseq_native breakseq_native [breakseq_native ...]

                        BreakSeq native GFF files (default: )

  --cnvnator_vcf cnvnator_vcf [cnvnator_vcf ...]

                        VCF file or dir for CNVnator VCFs (default: )

  --cnvnator_native File list [File list ...]

                        CNVnator native files (default: )

  --gatk_vcf file [file ...]

                        VCF file or dir for gatk VCFs (default: )

  --manta_vcf MANTA_VCF [MANTA_VCF ...]

                        VCF file or dir for Manta VCFs (default: )

  --lumpy_vcf LUMPY_VCF [LUMPY_VCF ...]

                        VCF file or dir for Lumpy VCFs (default: )

  --cnvkit_vcf CNVKIT_VCF [CNVKIT_VCF ...]

                        VCF file or dir for CNVkit VCFs (default: )

  --wham_vcf WHAM_VCF [WHAM_VCF ...]

                        VCF file or dir for WHAM VCFs (default: )

  --mean_read_length MEAN_READ_LENGTH

                        Mean read length (default: 100)

Reference options:

  --reference reference

                        Reference file (default: None)

  --chromosomes CHROMOSOMES [CHROMOSOMES ...]

                        Chromosome list to process. If unspecified, then all

                        chromosomes will be considered. (default: )

  --gaps gaps           Gap bed file (default: None)

  --filter_gaps         Filter out gaps (default: False)

  --keep_standard_contigs

                        Keep only the major contigs + MT (default: False)

Input BAM options:

  --bams BAMS [BAMS ...]

                        BAMs (default: [])

  --isize_mean ISIZE_MEAN

                        Insert size mean (default: 350.0)

  --isize_sd ISIZE_SD   Insert size standard deviation (default: 50.0)

Tool output merging options:

  --wiggle WIGGLE       Wiggle for interval overlap (default: 100)

  --inswiggle INSWIGGLE

                        Wiggle for insertions, overides wiggle (default: 100)

  --minsvlen MINSVLEN   Minimum length acceptable to be an SV (default: 50)

  --maxsvlen MAXSVLEN   Maximum length SV to report (default: 1000000)

  --overlap_ratio OVERLAP_RATIO

                        Reciprocal overlap ratio (default: 0.5)

Soft-clip detection options:

  --min_avg_base_qual MIN_AVG_BASE_QUAL

                        Minimum average base quality (default: 20)

  --min_mapq MIN_MAPQ   Minimum MAPQ (default: 5)

  --min_soft_clip MIN_SOFT_CLIP

                        Minimum soft-clip (default: 20)

  --max_nm MAX_NM       Maximum number of edits (default: 10)

  --min_matches MIN_MATCHES

                        Mininum number of matches (default: 50)

  --min_support_ins MIN_SUPPORT_INS

                        Minimum read support for calling insertions using

                        soft-clips (including neighbors) (default: 15)

  --min_support_frac_ins MIN_SUPPORT_FRAC_INS

                        Minimum fraction of reads supporting insertion using

                        soft-clips (including neighbors) (default: 0.05)

  --max_ins_intervals MAX_INS_INTERVALS

                        Maximum number of insertion intervals to generate

                        (default: 10000)

  --mean_read_coverage MEAN_READ_COVERAGE

                        Mean read coverage (default: 50)

  --min_ins_cov_frac MIN_INS_COV_FRAC

                        Minimum read coverage around the insertion breakpoint.

                        (default: 0.5)

  --max_ins_cov_frac MAX_INS_COV_FRAC

                        Maximum read coverage around the insertion breakpoint.

                        (default: 1.5)

  --sc_other_scale SC_OTHER_SCALE

                        Control degree of incorporation of breakpoints from

                        other methods. (default: 5)

Assembly options:

  --spades SPADES       Path to SPAdes executable (default: None)

  --spades_options SPADES_OPTIONS

                        Options for SPAdes (default: )

  --spades_timeout SPADES_TIMEOUT

                        Maximum time (in seconds) for running SPAdes on an

                        interval (default: 300)

  --disable_assembly    Disable assembly (default: False)

  --svs_to_assemble {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]

                        SVs to assemble (default: ['INS', 'INV', 'DUP'])

  --svs_to_softclip {DUP,INV,DEL,INS} [{DUP,INV,DEL,INS} ...]

                        SVs to soft-clip (default: ['INS', 'INV', 'DUP'])

  --extraction_max_read_pairs EXTRACTION_MAX_READ_PAIRS

                        Maximum number of pairs to extract for assembly

                        (default: 10000)

  --spades_max_interval_size SPADES_MAX_INTERVAL_SIZE

                        Maximum SV length for assembly (default: 50000)

  --assembly_max_tools ASSEMBLY_MAX_TOOLS

                        Skip assembly if more than this many tools support a

                        call (default 1) (default: 1)

  --assembly_pad ASSEMBLY_PAD

                        Padding base pairs to use for assembling breakpoint

                        with Spades and AGE (default: 500)

  --stop_spades_on_fail

                        Abort on SPAdes failure (default: False)

  --age AGE             Path to AGE executable (default: None)

  --age_timeout AGE_TIMEOUT

                        Maximum time (in seconds) for running AGE on an

                        assembled contig (default: 300)

  --min_inv_subalign_len MIN_INV_SUBALIGN_LEN

                        Minimum length of inversion sub-alginment (default:

                        50)

  --min_del_subalign_len MIN_DEL_SUBALIGN_LEN

                        Minimum length of deletion sub-alginment (default: 50)

  --age_window AGE_WINDOW

                        Window size for AGE to merge nearby breakpoints

                        (default: 20)

  --boost_sc            Use soft-clips for improving breakpoint detection

                        (default: False)

Genotyping options:

  --gt_window GT_WINDOW

                        Window for genotyping (default: 100)

  --gt_normal_frac GT_NORMAL_FRAC

                        Min. fraction of reads supporting reference for

                        genotyping (default: 0.05)

Output options:

  --svs_to_report {INV,CTX,INS,DEL,ITX,DUP} [{INV,CTX,INS,DEL,ITX,DUP} ...]

                        SV types to report (default: set(['INV', 'CTX', 'INS',

                        'DEL', 'ITX', 'DUP']))

  --enable_per_tool_output

                        Enable output of merged SVs for individual tools

                        (default: False)

Running environment options:

  --workdir WORKDIR     Scratch directory for working (default: work)

  --num_threads NUM_THREADS

                        Number of threads to use (default: 1)

  --outdir OUTDIR       Output directory (default: None)

Other options:

  --version             show program's version number and exit

——

テストラン

git clone https://github.com/bioinform/metasv.git
cd metasv/test/
./test_run.sh

 test_run.sh実行前

実行後。out/にvariants.vcfができている。

variants.vcfの中身。

動作している。

ラン

いくつかのSVコーラーは以前まとめて紹介してます（リンク）。ここでは各ツールの SVコールコマンドと、それを入力としてmetaSVを実行する手順を確認する。

１、mapping

ペアエンドのfastqからbamを作成する。リファレンスはhuman_g1k_v37_decoy.faとする。minimap2でマッピングし（*1）、samtoolsにパイプしてソートし、BAM出力（*2）（*6, *7）。

#normal fastq
minimap2 -ax sr -R "@RG	ID:X	LB:Y	SM:NA12877	PL:ILLUMINA" 
-t 16 human_g1k_v37_decoy.fa pair1.fg pair2.fq 
|samtools sort -@ 16 -m 4g -O BAM - > sample.bam
samtools index sample.bam

PCR duplicationを除き、リアライメントまで実行したいときは、こちらを参照。

２、Breakdancer-maxでSVコール（リンク）。dokcer imageを使う。

#コンテナ内で作業する
docker pull molecular/breakdancer #pull
docker run --rm -itv /Volumes/10TB/test/:/data/ molecular/breakdancer

> bam2cfg.pl sample.bam > breakdancer.cfg
> breakdancer-max breakdancer.cfg > breakdancer_output

３、PindelでSVコール（リンク）。dokcer imageを使う。

#コンテナ内で作業する
docker pull opengenomics/pindel
docker run -m "50g" --rm -itv /Volumes/10TB/test/:/data/ opengenomics/pindel


> #configファイル作成
> printf "sample.bam_350_sample1" > config
> sed -e 's/_/ /g' config > pindel_config && rm config #tabを入れたかったので少し遠回りした
> #pindel実行
> pindel -f human_g1k_v37_decoy.fa -T 10 -i pindel_config -o pindel
> cat pindel_SI pindel_D pindel_LI pindel_TD > pindel_merge
> pindel2vcf -R human_g1k_v37_decoy.fa -r human_g1k_v37_decoy.fa -p pindel_merge -v pindel.vcf -d 20180801 -e 4
> exit #コンテナを抜ける

vcffilter -f " SVTYPE = INS | SVTYPE = DEL " pindel.vcf > pindel.indel.vcf

４、MantaでSVコール（リンク）。dokcer imageを使う。

docker pull eitanbanks/manta-1.0.3
#configファイル作成
docker run --rm -itv /Volumes/10TB/test/:/data eitanbanks/manta-1.0.3 
/bin/manta/bin/configManta.py 
--normalBam=/data/sample.bam 
--referenceFasta=/data/human_g1k_v37_decoy.fa 
--runDir=/data/output

#configファイルを実行
docker run --rm -itv /Volumes/10TB/test/:/data eitanbanks/manta-1.0.3   
data/output/runWorkflow.py -m local -j 8 --memGb=20

５、LumpyでSVコール（リンク）。Lumpyはbrewで導入した。

#normal bam
samtools view -@ 12 -b -F 1294 sample.bam > sample.discordants.unsorted.bam
samtools view -@ 12 -h sample.bam | extractSplitReads_BwaMem -i stdin | samtools view -Sb - > sample.splitters.unsorted.bam
samtools sort -@ 12 sample.discordants.unsorted.bam > sample.discordants.bam
samtools sort -@ 12 sample.splitters.unsorted.bam > sample.splitters.bam

#lumpyを実行
lumpyexpress -B sample.bam -S sample.splitters.bam -D sample.discordants.bam -o lumpy.vcf

６、CnvkitでSVコール（リンク）。Cnvkitはcondaで導入した。dockerイメージ使うなら"docker pull etal/cnvkit"。

cnvkit.py batch sample.bam --normal \
 --targets my_baits.bed --fasta human_g1k_v37_decoy.fa \
 --output-reference sample.cnn --output-dir cnvkit/

#.cnsからvcfを出力
cnvkit.py export vcf cnvkit/sample.cns -i "SampleID" -o sample.cnv.vcf

７、CnvnatorでSVコール（リンク）。dokcer imageを使う。メモリ要求量が大きいようなので、メモリlimitは100gとした（*4）。

#オフィシャルらしきものもあるが、timothyjamesbeckerさんのをpullした（リンク）。
docker pull timothyjamesbecker/cnvnator
docker run -m "100g" --rm -itv /Volumes/10TB/test/:/data/ timothyjamesbecker/cnvnator
#作成中

SVコールが終わったら、MetaSV実行する。

方法１

Only merging output of 5 SV detectors without further sof-clips based analysis or local assembly.

SV callerの結果のマージだけ実行。soft clipやローカルアセンブリを介した解析は行わない。

run_metasv.py --reference human_g1k_v37_decoy.fa 
--lumpy_vcf lumpy_tumor_normal.vcf 
--cnvkit_vcf Sample.cnv.vcf 
--breakdancer_native breakdancer_output 
--manta_vcf manta.vcf
--pindel_vcf pindel.indel.vcf
--outdir metasv_out --sample HG005 --disable_assembly --filter_gaps --keep_standard_contigs --num_threads 12 --workdir work

出力ディレクトリ。統合されフィルタリングされたVCFファイルが出力される。

方法２

Complete run of MetaSV using all 4 SV detectors, soft-clips based analysis to enhance insertion detection, and local assembly to improve breakpoint resolution.

全解析。

run_metasv.py --reference human_g1k_v37_decoy.fa 
--lumpy_vcf lumpy_tumor_normal.vcf 
--cnvkit_vcf Sample.cnv.vcf 
--breakdancer_native breakdancer_output 
--manta_vcf manta.vcf
--pindel_vcf pindel.indel.vcf
--spades SPAdes/spades.py  --age AGE/age_align 
--min_support_ins 2 --max_ins_intervals 500000 --isize_mean 500 
--isize_sd 150 --num_threads 12 --workdir work --outdir out  
--sample HG005 --bam alignments.bam 

mac osでは、アセンブリしたcontigのindex作成あたりでエラーを起こした。やはりlinux機でテストした方が良かったかもしれない。

他の方法は著者らが準備した説明ページで確認して下さい。BreakSeq2のコマンドはわからなかったので載せてません。GATKについては、バージョンが異なると内容が変わり、余計な混乱を招くので載せてません（GATK３と最新のGATK４の違いとか）。

MetaSV 

より新しい手法も発表されてきています。また紹介します。

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

引用

MetaSV: an accurate and integrative structural-variant caller for next generation sequencing
Marghoob Mohiyuddin, John C. Mu, Jian Li, Narges Bani Asadi, Mark B. Gerstein,Alexej Abyzov, Wing H. Wong, and Hugo Y.K. Lam

Bioinformatics. 2015 Aug 15; 31(16): 2741–2744. 

*1 elprepでの簡単フィルタイングを考えたが、elprepの一部タスクは全データをメモリ内に収納して動作するためメモリ要求量が高く、chrごとにsplitしてもメモリが足りないのでやめた（ヒト50xWGSでも、おそらくは384GBくらいメモリが確保されていれば、splitで動く）

elprep split

　　 ↓

elprep filter --mark-duplicates --remove-duplicates --filter-mapping-quality 1 --clean-sam

*2 mac osでランしたかったので、処理速度の低下には目をつぶり、dokcer imagesを使っている。

*3 使用スレッド数を増やす場合、利用可能なメモリー量などに応じて慎重に行って下さい。

*4 実メモリに応じて変えて下さい。

*5 新しくても動くかもしれない。未検討。

*6 メモリ追加。指定値xスレッド数、だけメモリを食います。

*7 sort時にtoo many  filesのエラーが出たら、このページを参考にopen filesの上限を上げて下さい（現在の値を確認するには"ulimit -n"）。

コピー数変化の検出と可視化ツール CNVkit

2020 3/22 実行例追記

2022/04/29 インストール手順修正

　コピー数変化は、ガンを含む多くの疾患の有用な診断指標である。ゲノム全体のコピー数解析のゴールドスタンダードは、 array comparative genomic hybridization（array CGH）である[論文より　ref.1,2]。より最近では、全ゲノムシーケンシングデータからコピー数の情報を得る方法が開発されている（非特許文献4参照）。臨床用途では、エキソンや疾患関連遺伝子のセットなどのシーケンス解析では、興味のある領域をより高いカバレッジでシーケンスして、バリアントコールの感受性を高めることが好まれる[ref.5]。これらのデータを解析するため、 CNVer [ref.6]、ExomeCNV [ref.7]、exomeCopy [ref.8]、CONTRA [ref.9]、CoNIFER [ref.10]、ExomeDepth [ref.11]、VarScan 2 [ref.12]、XHMM [ref.13]、ngCGH [ref.14]、EXCAVATOR [ref.15]、CANO​​ES [ref.16]、PatternCNV [ref.17]、CODEX [ref.18]、Control-FREEC [ref.19]およびcn.MOPS [ref.20]などのコピー数を解析するツールも開発されてきた。しかしながら、これらのアプローチは、遺伝子間および通常はイントロン領域のシーケンスリードを使用せず、コピー数を推測する可能性を制限する。

　ターゲットエンリッチメントでは、標的領域はハイブリダイゼーションによりキャプチャされる。しかし、かなりの量のオフターゲットDNAがライブラリーに残っており、このDNAはシーケンスされてリードのかなりの部分を占めている。したがって、オフターゲットのリードは、全ゲノムに渡り非常に低カバレッジのシーケンスを提供する。ターゲット外のリードだけではSNVおよび他の小さなバリアントをコールする十分なカバレッジは得られないが、最近cnvOffSeq [ref.21]およびCopywriteR [ref.22]によって実証されたように、ラージスケールではコピー数推測に有用な情報を提供する。

　著者らは、ターゲットシーケンスでコピー数変化および改変の分析のための計算方法を開発した。このツールキットはCNVkitと呼ばれ、CNV検出のパイプラインを実装している。オンとオフの両方のターゲットシーケンシングリードを利用し、一連の修正を適用してコピー数のコール精度を向上させている。オンとオフのターゲット領域のビニングされたリードデプスを比較し、異なる解像度であるにもかかわらず、コピー数の類似した推定値を提供することを見出した。真のコピー数変化によって駆動されていないであろうビニングされたリードカウントの変動を補正するアルゴリズムをいくつか評価した。最後に、CNVkitメソッドと2つの競合するCNVコーラーの結果をarray CGHと比較し、CNVkitがarray CGHと最もよく一致することを確認した。

マニュアル（丁寧にまとめられています。ボリュームがあります。）

https://cnvkit.readthedocs.io/en/stable/

Biostars quwestion

https://www.biostars.org/t/CNVkit/

インストール

mac os10.13のanaconda2-4.3.0でテストした。

依存

• python >=3.5

Github

#bioconda (コメントいただきました通り修正しました。)
mamba create -n cnvkit -c conda-forge -c bioconda python=3.9 cnvkit -y
conda activate cnvkit

pipでも導入できる。Dockerコンテナも準備されている。
#dockerイメージ
docker pull etal/cnvkit

> cnvkit.py -h

]$ cnvkit.py -h

usage: cnvkit.py [-h]

                 {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,import-picard,import-seg,import-theta,import-rna,export,version}

                 ...

CNVkit, a command-line toolkit for copy number analysis.

positional arguments:

  {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,import-picard,import-seg,import-theta,import-rna,export,version}

                        Sub-commands (use with -h for more info)

    batch               Run the complete CNVkit pipeline on one or more BAM

                        files.

    target              Transform bait intervals into targets more suitable

                        for CNVkit.

    access              List the locations of accessible sequence regions in a

                        FASTA file.

    antitarget          Derive off-target ("antitarget") bins from target

                        regions.

    autobin             Quickly calculate reasonable bin sizes from BAM read

                        counts.

    coverage            Calculate coverage in the given regions from BAM read

                        depths.

    reference           Compile a coverage reference from the given files

                        (normal samples).

    fix                 Combine target and antitarget coverages and correct

                        for biases. Adjust raw coverage data according to the

                        given reference, correct potential biases and re-

                        center.

    segment             Infer copy number segments from the given coverage

                        table.

    call                Call copy number variants from segmented log2 ratios.

    diagram             Draw copy number (log2 coverages, CBS calls) on

                        chromosomes as a diagram. If both the raw probes and

                        segments are given, show them side-by-side on each

                        chromosome (segments on the left side, probes on the

                        right side).

    scatter             Plot probe log2 coverages and segmentation calls

                        together.

    heatmap             Plot copy number for multiple samples as a heatmap.

    breaks              List the targeted genes in which a copy number

                        breakpoint occurs.

    genemetrics         Identify targeted genes with copy number gain or loss.

    sex                 Guess samples' sex from the relative coverage of

                        chromosomes X and Y.

    metrics             Compute coverage deviations and other metrics for

                        self-evaluation.

    segmetrics          Compute segment-level metrics from bin-level log2

                        ratios.

    import-picard       Convert Picard CalculateHsMetrics tabular output to

                        CNVkit .cnn files. The input file is generated by the

                        PER_TARGET_COVERAGE option in the CalculateHsMetrics

                        script in Picard tools. If 'antitarget' is in the

                        input filename, the generated output filename will

                        have the suffix '.antitargetcoverage.cnn', otherwise

                        '.targetcoverage.cnn'.

    import-seg          Convert a SEG file to CNVkit .cns files.

    import-theta        Convert THetA output to a BED-like, CNVkit-like

                        tabular format. Equivalently, use the THetA results

                        file to convert CNVkit .cns segments to integer copy

                        number calls.

    import-rna          Convert a cohort of per-gene log2 ratios to CNVkit

                        .cnr format.

    export              Convert CNVkit output files to another format.

    version             Display this program's version.

optional arguments:

  -h, --help            show this help message and exit

See the online manual for details: https://cnvkit.readthedocs.io

たくさんのコマンドがある。

例えばcnvkitの全パイプラインを実行するbatchコマンドのヘルプを表示する。

> cnvkit.py batch -h

 $ cnvkit.py batch -h

usage: cnvkit.py batch [-h] [-m {hybrid,amplicon,wgs}] [-y] [-c]

                       [--drop-low-coverage] [-p [PROCESSES]]

                       [--rlibpath DIRECTORY] [-n [FILES [FILES ...]]]

                       [-f FILENAME] [-t FILENAME] [-a FILENAME]

                       [--annotate FILENAME] [--short-names]

                       [--target-avg-size TARGET_AVG_SIZE] [-g FILENAME]

                       [--antitarget-avg-size ANTITARGET_AVG_SIZE]

                       [--antitarget-min-size ANTITARGET_MIN_SIZE]

                       [--output-reference FILENAME] [-r REFERENCE]

                       [-d DIRECTORY] [--scatter] [--diagram]

                       [bam_files [bam_files ...]]

positional arguments:

  bam_files             Mapped sequence reads (.bam)

optional arguments:

  -h, --help            show this help message and exit

  -m {hybrid,amplicon,wgs}, --method {hybrid,amplicon,wgs}

                        Sequencing protocol: hybridization capture ('hybrid'),

                        targeted amplicon sequencing ('amplicon'), or whole

                        genome sequencing ('wgs'). Determines whether and how

                        to use antitarget bins. [Default: hybrid]

  -y, --male-reference, --haploid-x-reference

                        Use or assume a male reference (i.e. female samples

                        will have +1 log-CNR of chrX; otherwise male samples

                        would have -1 chrX).

  -c, --count-reads     Get read depths by counting read midpoints within each

                        bin. (An alternative algorithm).

  --drop-low-coverage   Drop very-low-coverage bins before segmentation to

                        avoid false-positive deletions in poor-quality tumor

                        samples.

  -p [PROCESSES], --processes [PROCESSES]

                        Number of subprocesses used to running each of the BAM

                        files in parallel. Without an argument, use the

                        maximum number of available CPUs. [Default: process

                        each BAM in serial]

  --rlibpath DIRECTORY  Path to an alternative site-library to use for R

                        packages.

To construct a new copy number reference:

  -n [FILES [FILES ...]], --normal [FILES [FILES ...]]

                        Normal samples (.bam) used to construct the pooled,

                        paired, or flat reference. If this option is used but

                        no filenames are given, a "flat" reference will be

                        built. Otherwise, all filenames following this option

                        will be used.

  -f FILENAME, --fasta FILENAME

                        Reference genome, FASTA format (e.g. UCSC hg19.fa)

  -t FILENAME, --targets FILENAME

                        Target intervals (.bed or .list)

  -a FILENAME, --antitargets FILENAME

                        Antitarget intervals (.bed or .list)

  --annotate FILENAME   Use gene models from this file to assign names to the

                        target regions. Format: UCSC refFlat.txt or

                        ensFlat.txt file (preferred), or BED, interval list,

                        GFF, or similar.

  --short-names         Reduce multi-accession bait labels to be short and

                        consistent.

  --target-avg-size TARGET_AVG_SIZE

                        Average size of split target bins (results are

                        approximate).

  -g FILENAME, --access FILENAME

                        Regions of accessible sequence on chromosomes (.bed),

                        as output by the 'access' command.

  --antitarget-avg-size ANTITARGET_AVG_SIZE

                        Average size of antitarget bins (results are

                        approximate).

  --antitarget-min-size ANTITARGET_MIN_SIZE

                        Minimum size of antitarget bins (smaller regions are

                        dropped).

  --output-reference FILENAME

                        Output filename/path for the new reference file being

                        created. (If given, ignores the -o/--output-dir option

                        and will write the file to the given path. Otherwise,

                        "reference.cnn" will be created in the current

                        directory or specified output directory.)

To reuse an existing reference:

  -r REFERENCE, --reference REFERENCE

                        Copy number reference file (.cnn).

Output options:

  -d DIRECTORY, --output-dir DIRECTORY

                        Output directory.

  --scatter             Create a whole-genome copy ratio profile as a PDF

                        scatter plot.

  --diagram             Create an ideogram of copy ratios on chromosomes as a

                        PDF.

——

試せてませんが、galaxyやDNA Nexusでも利用できるようです。

Galaxy | (sandbox for testing) | Tool Shed

https://platform.dnanexus.com/login

実行方法

--CNVの検出--

様々なサブコマンドがあるが、ここでは全パイプラインを動かすbatchコマンドと、出力を変換するexportコマンドを紹介する（Quick startより）。 

batch: Run the complete CNVkit pipeline on one or more BAM files.

パターンA:　turmor-normal

cnvkit.py batch *Tumor.bam --normal *Normal.bam \
 --targets my_baits.bed --fasta hg19.fasta \
 --access data/access-5kb-mappable.hg19.bed \
 --output-reference my_reference.cnn --output-dir example/

パターンB:　比較対象のbamがない場合、baitリファレンスを指定することでランできる。その場合、"--normal"フラグを引数なしで立てる。

cnvkit.py batch *Tumor.bam --normal -t my_baits.bed -f hg19.fasta \
 --access data/access-5kb-mappable.hg19.bed \
 --output-reference my_flat_reference.cnn -d example/

 ランが終わると、出力ディレクトリにいくつかのファイルができる。

• Target and antitarget bin-level coverages (.cnn)

• Copy number reference profile (.cnn)

• Bin-level log2 ratios (.cnr)

• Segmented log2 ratios (.cns)

詳細はフォーマット解説参照。

https://cnvkit.readthedocs.io/en/stable/fileformats.html

 出力をvcfに変換する。

export: Convert CNVkit output files to another format.

cnvkit.py export vcf example/sample.cns -i "SampleID" -o sample.cnv.vcf

 --CNVの可視化--

Bin-level log2 ratios (.cnr)ファイルとsegmentation callsをプロットする（おそらくhumanゲノムを想定しているので注意）

scatter（リンク）: Plot probe log2 coverages and segmentation calls
together.

cnvkit.py scatter example/sample.cnr -s example/sample.cns

diagram（リンク）: Draw copy number (log2 coverages, CBS calls) on chromosomes as a diagram.

cnvkit.py diagram -s example/sample.cns example/sample.cnr

heatmap（リンク）: Plot copy number for multiple samples as a heatmap.

cnvkit.py heatmap example/*.cns

heatmap機能は縦にサンプルを並べて描画するので（リンク）、familyやpopulation解析など、一定の数のサンプルを比較する時に使えるんじゃないかと思います。

追記

コンティグごとのカバレッジを出力する。

#contig.faのindexからbedを作る。
samtools faidx contig.fa
awk '{print $1 "\t0\t" $2}' contig.fa.fai > contig.bed

#cnvkit.pyのcoverageコマンド実行
cnvkit.py coverage -o out input.bam contig.bed

出力

引用

CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing

Eric Talevich, A. Hunter Shain, Thomas Botton, Boris C. Bastian

PLoS Comput Biol. 2016 Apr; 12(4): e1004873.

オリジナルfastqと比較してbamのリード情報が完全に同じかどうか調べる BamHash

　（ゲノム）リシーケンシングプロジェクトは、既知ゲノムを有する種の個体のシーケンシング解析であり、大量のraw シーケンシングリードを生成し、その後、これらはリファレンスゲノムにアライメントされる。シーケンシングコストが減少し、現在のシーケンシング技術のスループットが増加し続けているため、データストレージが問題になっている。

　Rawシーケンシングリードは、一般的に圧縮されたFASTQファイル形式で保存される。マッピング後に得られたアライメント情報はBAMファイルに格納される（Li et al、2009）。このBAMファイルはソートされ、さらに処理されるが、最も重要なのは、FASTQファイルのオリジナル情報が全て含まれていることである。ソートされたBAMファイルは、ソートされていないBAMファイルと比較してより優れた圧縮を提供し、ゲノム領域に対するランダムルックアップを可能にする。この理由のために、ほとんど全てのポストアラインメント解析が可能である。例えばバリアントコール、リアライメント、ローカルアセンブリは、元のFASTQファイルではなく、ソートされたBAMファイルを使って実行される。

　BAMファイルにはFASTQファイルのすべての情報が含まれているため、アライメント後にFASTQファイルを削除することには正当性がある。つまりFASTQファイルの内容はBAMファイルから再生成することができる。

　しかし、FASTQファイルを削除する前に、データの損失がないこと、すなわちFASTQファイルのシーケンスがBAMファイルのシーケンスとまったく同じであることを確認する必要がある。 2つのファイルは様々な理由により異なる場合がある。アライメントパイプラインにエラーがあると、一貫性のないファイルが生成される可能性がある。アライメントパイプラインは十分にテストされたツールに基づいているが、通常の状態で動作することを意図しており、ハードウェアの故障やディスク領域がなくなった場合に動作が予測できない可能性がある。したがって、パイプライン全体の出力を独立して検証できることが重要になる。

　（この論文では）FASTQとBAMファイル間のデータの整合性を検証するためのツールであるBamHashを紹介する。このプログラムは両方のfastqの配列とリード名から64ビットのfingerprint（用語）を計算する。この方法は入力の変化に対して非常に敏感であるため、１ヌクレオチドの変化でも異なるfingerprintを生じる。偶然同じfingerprintを生成する確率は天文学的に小さい。このツールの役割は、異なるフィンガープリントを持つFASTQファイルとBAMファイルにフラグを立てて、FASTQファイルの削除が安全でないとマークをつけることである。

　BamHashは、ファイルのfingerprintを計算するmd5sumプログラムと同じ役割を果たす。フォーマットと順序が異なるため、FASTQとBAMファイルのmd5sumのfingerprint（Rivest、1992）を比較しても同等の結果は得られない。著者らの方法は高速かつメモリ効率的である。30xヒトゲノムシーケンシングのBAMファイルのfingerprintを30分で計算することができる。

BamHashに関するツイート。

インストール

ubuntu18.04でテストした。

 Github

git clone https://github.com/DecodeGenetics/BamHash.git
cd BamHash/
make

> ./bamhash_checksum_bam -h

$ ./bamhash_checksum_bam -h

bamhash_checksum_bam - Checksum of a bam file

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

SYNOPSIS

    bamhash_checksum_bam [OPTIONS] <in.bam> <in2.bam> ...

DESCRIPTION

    Program for checksum of sequence reads.

    -h, --help

          Displays this help message.

    --version

          Display version information

  Options:

    -d, --debug

          Debug mode. Prints full hex for each read to stdout

    -R, --no-readnames

          Do not use read names as part of checksum

    -Q, --no-quality

          Do not use read quality as part of checksum

    -P, --no-paired

          Bam files were not generated with paired-end reads

VERSION

    bamhash_checksum_bam version: 1.1

    Last update May 2015

—

> ./bamhash_checksum_fasta -h

$ ./bamhash_checksum_fasta -h

bamhash_checksum_fasta - Checksum of a set of fasta files

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

SYNOPSIS

    bamhash_checksum_fasta [OPTIONS] <in1.fasta> [in2.fasta ... ]

DESCRIPTION

    Program for checksum of sequence reads.

    -h, --help

          Displays this help message.

    --version

          Display version information

  Options:

    -d, --debug

          Debug mode. Prints full hex for each read to stdout

    -R, --no-readnames

          Do not use read names as part of checksum

VERSION

    bamhash_checksum_fasta version: 1.1

    Last update May 2015

——

> ./bamhash_checksum_fastq -h

$ ./bamhash_checksum_fastq -h

bamhash_checksum_fastq - Checksum of a set of fastq files

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

SYNOPSIS

    bamhash_checksum_fastq [OPTIONS] <in1.fastq.gz> [in2.fastq.gz ... ]

DESCRIPTION

    Program for checksum of sequence reads.

    -h, --help

          Displays this help message.

    --version

          Display version information

  Options:

    -d, --debug

          Debug mode. Prints full hex for each read to stdout

    -R, --no-readnames

          Do not use read names as part of checksum

    -Q, --no-quality

          Do not use read quality as part of checksum

    -P, --no-paired

          List of fastq files are not paired-end reads

VERSION

    bamhash_checksum_fastq version: 1.1

    Last update May 2015

kazu@ubuntu:~/BamHash$

——

ラン

ペアエンドのbamファイル（シングルエンドは"--no-paired"をつける）

bamhash_checksum_bam [OPTIONS] <in.bam> <in2.bam> ...

ペアエンドのfastq

bamhash_checksum_fastq [OPTIONS] <in1.fastq.gz> [in2.fastq.gz ... ]

ペアエンドのFASTAファイル（シングルエンドは"--no-paired"をつける）

bamhash_checksum_fastq [OPTIONS] <in1.fastq.gz> [in2.fastq.gz ... ]

ペアエンドのFASTAからbamを作成した時は、bamのchecksumを"--no-quality"フラグをつけて行う。

テスト

とあるbamファイルのハッシュ値（fastq部分だけのハッシュ値）を調べる。

$ bamhash_checksum_bam recal_reads.bam 

c7d51310b2a044d3 508742

次にそのfastqのchecksumを調べる。

$ bamhash_checksum_fastq R1.fastq R2.fastq 

1e5b86c9181ca09f 254371

異なると判定された。これは、sam=>bamにするときにelprep（リンク）で不要なリード情報を除いたbamを使ったためである。次は、先ほどと同じbamだがフィルタリングせず作ったbamのハッシュ値。

$ bamhash_checksum_bam R1R2_sorted.bam 

1e5b86c9181ca09f 508742

合致している。

引用
BamHash: a checksum program for verifying the integrity of sequence data
Óskarsdóttir A, Másson G, Melsted P

Bioinformatics. 2016 Jan 1;32(1):140-1