macでインフォマティクス

macでインフォマティクス

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

指定されたカバレッジになるようにリードをランダムサブサンプリングする rasusa

 

レポジトリより

 私(オーサー)の要件を満たすリードのサブサンプリング用ツールは見つかりませんでした。私が見つけたツールは、サブサンプルするリードの数や割合を指定するだけのものや、カバレッジに合わせてサブサンプルする場合、すべてのリードが同じサイズであると仮定するもの(Illumina など)など、どれも不十分でした。私は主にロングリードのデータを扱っているので、ファイルを特定のカバレッジにサブサンプルしたい場合、リードの長さが考慮されていないため、これが問題となっていました。

 私がしばらく使っていた回避策は、filtlong(紹介);サンプルの理論的なカバレッジを達成するために必要な塩基数を計算するだけの簡単なもの、を使うことでした。例えば、大腸菌のサンプルで500万リードのfastqがあり、それを50倍のカバレッジになるようにサブセットしたいとします。サンプルのゲノムの予想サイズである460万塩基対に希望するカバレッジを掛けるだけで、目標とする塩基数である2億3000万塩基対が得られます。filtlongでは、次のようにします。

ターゲット=230000000
$ filtlong --target_bases "$target" reads.fq > reads.50x.fq
しかし、これは技術的にはfiltlongの本来の機能ではなく、質の高いフィルタリングツールです。最終的に得られるのは、(理論上の)50倍のカバレッジで「最もスコアの高い」リードのサブセットです。あなたの状況によっては、これはあなたが望むものかもしれません。しかし、データセットの中で最も良い/長いリードに偏ってしまうので、データセット全体を公平に表すことはできません。また、より長く、より質の高いリードを生成するゲノムの領域が優先される可能性もあります。De Maioらは、ナノポアのリードをランダムにサブサンプリングすることで、フィルタリングした場合よりも優れたゲノムアセンブリが得られることを発見しました。

つまり、あなたの状況によっては、リードの偏りのないサブサンプルが必要になるかもしれないのです。そのような場合、rasusaはあなたをサポートします。

 

インストール

Github

#conda (link)
conda install -c bioconda -y rasusa

#cargo 
cargo install rasusa

#brew
brew tap brewsci/bio
brew install rasusa

#docker
docker pull quay.io/mbhall88/rasusa

> rasusa -h

rasusa 0.6.0

Randomly subsample reads to a specified coverage.

 

USAGE:

    rasusa [FLAGS] [OPTIONS] --bases <bases> --coverage <FLOAT> --genome-size <size|faidx> --input <input>...

 

FLAGS:

    -h, --help       Prints help information

    -V, --version    Prints version information

    -v               Switch on verbosity.

 

OPTIONS:

    -b, --bases <bases>               Explicitly set the number of bases required e.g., 4.3kb, 7Tb, 9000, 4.1MB

    -l, --compress-level <1-9>        Compression level to use if compressing output [default: 6]

    -c, --coverage <FLOAT>            The desired coverage to sub-sample the reads to

    -g, --genome-size <size|faidx>    Genome size to calculate coverage with respect to. e.g., 4.3kb, 7Tb, 9000, 4.1MB

    -i, --input <input>...            The fast{a,q} file(s) to subsample

    -o, --output <output>...          Output filepath(s); stdout if not present

    -O, --output-type <u|b|g|l>       u: uncompressed; b: Bzip2; g: Gzip; l: Lzma

    -s, --seed <INT>                  Random seed to use.

 

 

実行方法

リードとゲノムサイズ、カバレッジ(リードデプス)を指定する。

rasusa -i in.fq.gz -c 30 -g 4.6mb -o out.fq

#gzipped fastqとして保存
rasusa -i in.fq.gz -c 30 -g 4.6mb -O g -o out.fq.gz
  • -g   Genome size to calculate coverage with respect to. e.g., 4.3kb, 7Tb, 9000, 4.1MB.
  • -i    The fast{a,q} file(s) to subsample
  • -o   Output filepath(s); stdout if not present
  • -O <u|b|g|l>       u: uncompressed; b: Bzip2; g: Gzip; l: Lzma 
  • -l     Compression level to use if compressing output [default: 6]
  • -s    Random seed to use.

 

引用

Hall, Michael B. Rasusa: Randomly subsample sequencing reads to a specified coverage. (2019). doi:10.5281/zenodo.3731394

 

関連