k-merの起源となる配列を見つける Back to sequences

2024/10/09追記

　生のシーケンスデータの処理に特化したバイオインフォマティクスツールの大部分は、k-mersの概念を多用している。これにより、データの冗長性（ひいてはメモリの圧迫）を減らし、シーケンスエラーを破棄し、操作可能で容易に比較できる固定サイズのオブジェクトを処分することができる。欠点は、各k-merとそれが属する元の配列セットとの間のリンクが一般に失われることである。この文脈で考慮されるデータ量を考えると、この関連性を探し出すにはコストがかかる。この研究では、「back_to_sequences 」を紹介する。「back_to_sequences 」は、興味のあるk-merの集合にインデックスを付け、配列の集合をストリーミングし、インデックス付けされたk-merを少なくとも1つ含むものを抽出するように設計されたシンプルなツールである。さらに、配列中のk-merの出現数も提供する。back_to_sequencesは1ミリ秒あたり約200のショートリードをストリームし、数億のリード中のk-merを数分で検索できる。

Documentation

https://b2s-doc.readthedocs.io/en/latest/index.html

My first paper on @JOSS_TheOJ ! https://t.co/kbhzvWh0Xp
— Pierre Marijon 🏳️‍🌈 (@pierre_marijon) 2024年9月24日

Hey kmers' fans!
The "back_to_sequences" tool https://t.co/IgCdQbOOzd is now published in JOSS https://t.co/12CulElV0v.
Find sequences (reads, unitigs, genes) related to a set of kmers in large datasets, in a matter of seconds.
More to read in this thread: https://t.co/Eh7QI33y3w pic.twitter.com/OhK19xJEcF
— Pierre Peterlongo (@p_peterlongo) 2024年9月23日

インストール

ubuntu22.04LTSでテストした（rustc --version => rustc 1.78.0）。

To compile from source

Rust

Github

git clone https://github.com/pierrepeterlongo/back_to_sequences.git
cd back_to_sequences
RUSTFLAGS="-C target-cpu=native" cargo install --path .

> back_to_sequences -h

Back to sequences: find the origin of kmers

Usage: back_to_sequences [OPTIONS] --in-kmers <IN_KMERS>

Options:

--in-kmers <IN_KMERS>

Input fasta file containing the original kmers

Note: back_to_sequences considers the content as a set of kmers

This means that a kmer is considered only once,

even if it occurs multiple times in the file.

If the stranded option is not used (default), a kmer

and its reverse complement are considered as the same kmer.

--in-sequences <IN_SEQUENCES>

Input fasta or fastq [.gz] file containing the original sequences (eg. reads).

The stdin is used if not provided

(and if `--in_filelist` is not provided neither) [default: ]

--in-filelist <IN_FILELIST>

Input txt file containing in each line a path to a fasta or fastq [.gz] file

containing the original sequences (eg. reads).

Note1: if this option is used, the `--out_filelist` option must be used.

The number of lines in out_filelist must be the same as in_filelist

Note2: Incompatible with `--in_sequences` [default: ]

--out-sequences <OUT_SEQUENCES>

Output file containing the filtered original sequences (eg. reads).

It will be automatically in fasta or fastq format depending on the input file.

If not provided, only the in_kmers with their count is output [default: ]

--out-filelist <OUT_FILELIST>

Output txt file containing in each line a path to a fasta or fastq [.gz] file

that will contain the related output file from the input files list [default: ]

--out-kmers <OUT_KMERS>

If provided, output a text file containing the kmers that occur in the reads

with their

* number of occurrences

* their occurrence positions if the --output_kmer_positions option is used

Note: if `--in_filelist` is used the output counted kmers are

those occurring the last input file of that list [default: ]

--counted-kmer-threshold <COUNTED_KMER_THRESHOLD>

If out_kmers is provided, output only reference kmers whose number of occurrences

is at least equal to this value.

If out_kmers is not provided, this option is ignored [default: 0]

--output-kmer-positions

If out_kmers is provided, either only count their number of occurrences (default)

or output their occurrence positions (read_id, position, strand) if this option is used

--output-mapping-positions

If provided, output matching positions on sequences in the

out_sequence file(s)

-k, --kmer-size <KMER_SIZE>

Size of the kmers to index and search [default: 31]

-m, --min-threshold <MIN_THRESHOLD>

Output sequences are those whose ratio of indexed kmers is in ]min_threshold; max_threshold]

Minimal threshold of the ratio (%) of kmers that must be found in a sequence to keep it (default 0%).

Thus by default, if no kmer is found in a sequence, it is not output. [default: 0]

--max-threshold <MAX_THRESHOLD>

Output sequences are those whose ratio of indexed kmers is in ]min_threshold; max_threshold]

Maximal threshold of the ratio (%) of kmers that must be found in a sequence to keep it (default 100%).

Thus by default, there is no limitation on the maximal number of kmers found in a sequence. [default: 100]

--stranded

Used original kmer strand (else canonical kmers are considered)

--query-reverse

Query the reverse complement of reads. Useless without the --stranded option

--no-low-complexity

Do not index low complexity kmers (ie. with a Shannon entropy < 1.0)

-t, --threads <THREADS>

Number of threads

Note: if not provided, the number of threads is set to the number of logical cores [default: 0]

-h, --help

Print help

-V, --version

Print version

実行方法

探索するk-mer配列のfastaと探索対象のリード、出力として返すリード名、さらに任意でk-merカウント数を示したテキストを指定する。

back_to_sequences --in-kmers compacted_kmers.fasta --in-sequences reads.fasta --out-sequences filtered_reads.fasta  --out-kmers counted_kmers.txt

--in-kmers Input fasta file containing the original kmers
--in-sequences Input fasta or fastq [.gz] file containing the original sequences (eg. reads)
--out-sequences Output file containing the filtered original sequences (eg. reads).
--out-kmers If provided, output a text file containing the kmers that occur in the reads with their * number of occurrences or * their occurrence positions if the --output_kmer_positions option is used
-k Size of the kmers to index and search [default: 31]
-t Number of threads

出力例

filtered_reads.fastaは、入力にfastqフォーマットのリードを指定した時はfastqフォーマットとして出力される。

> head counted_kmers.txt

$ head counted_kmers.txt

CGCTCTATCTGAAACGCGGCGATACGATTTA 1

CGGAATTTGGAGAGAGTATTTTGGTCACGGC 1

ACCGACACGAAGCGACCCTCGCCCTGATCAT 2

CGCAACGTCGCGCGGATGACCGTGCTGCTCC 1

CTGCTCGATCTGTTCTTTGGTGACGGGCTTA 1

CCGCGCCAGCTTGATCGCCATAGTGCCAATC 1

ACATCGCCCACTTCGCGCAGACGAAAGCCCG 1

AACCTTCCGCAAAGGCAAGCGCCATACCGAC 1

ATTTTTCGGTCGCCATTGGCGCAAGCGCTGA 1

CGGCGGTACGTTCACATTCTCAGCGAGCTCC 1

compacted_kmers.fastaに50%以上、多くても70%以下のkmerが含まれる配列のみ出力するフィルターを指定。

back_to_sequences --in-kmers compacted_kmers.fasta --in-sequences reads.fasta --out-sequences filtered_reads.fasta  --out-kmers counted_kmers.txt --min-threshold 50 --max-threshold 70

--min-threshold Output sequences are those whose ratio of indexed kmers is in ]min_threshold
--max-threshold Output sequences are those whose ratio of indexed kmers is in ]min_threshold; max_threshold]

ほかにもstrand指定や逆相補の配列をターゲットにするオプションなどがある（Document参照）。

論文より

配列の集合から目的の一意なk-merを見つけるには、古典的なgrepや、The Platinum Searcher（PT）や The Silver Searcherのような最近のパターンマッチングツールを使うことができる。テストとして、MacBook、Apple M2 proを使用し、grep、The Platinum Searcher、The Silver Searcherを使用して、Q100Mデータセット（セクション2.1参照）の1つのk-merを照会した。grepは44秒を要し、ptは15秒、agは400秒後に停止した（result）。従って、単純に外挿すると、grepの場合は1台のコンピューターで100万個のk-merを検索するのに約500日かかり、back_to_sequencesを使った場合は5～6分である。