RNAseqのロングリードのアライメントの評価ツール RNAseqEval

ロングリードを使ったRNA seqはまだ情報が少ない。Evaluation of tools for long read RNA-seq splice-aware alignment.というタイトルのこの論文では、PacBioとONT Minionを使い、エラーの多いロングリードがアライナーの種類によってどう扱われれるのか、調べられた。その評価に使われたツールがGithubで公開されているので紹介する。

論文で使われたデータ

https://github.com/kkrizanovic/RNAseqEval/blob/master/RNAseq_benchmark/RNAseq_benchmark.md

インストール

Prerequisites

Python 2

本体　Github

git clone --recursive https://github.com/kkrizanovic/RNAseqEval.git

実行方法

リファレンスからシミュレートのRNA seqデータを発生(example_dataset/を使う）。

python generate_transcriptome.py example_dataset/dmelanogaster_chr4.gtf example_dataset/dmelanogaster_chr4_genome.fa transcriptome.fasta

gtfで指定した領域からロングリードがジェネレートされる。

> seqkit stats transcriptome.fasta

user$ seqkit stats transcriptome.fasta

file format type num_seqs sum_len min_len avg_len max_len

/Users/user/local/RNAseqEval/transcriptome.fasta FASTA DNA 334 1,336,895 69 4,002.7 27,666

マッピングを評価する。（リアルデータ向けだが、ここでは上で作ったシミュレートのrングリードを使う）

#最初にBBmap2でロングリードをアライメントする(条件は論文の 2.3 Evaluated RNA-seq toolsの通り)
bbmap.sh in=transcriptome.fasta out=mapped.sam ref=example_dataset/dmelanogaster_chr4_genome.fa fastareadlen=500

mapped.samができる。

mapping率を調べておく。

samtols view -bS mapped.sam > mapped.bam 
bamtools stats -in mapped.bam

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

Stats for BAM file(s):

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

Total reads: 2833

Mapped reads: 2833 (100%)

Forward strand: 1392 (49.1352%)

Reverse strand: 1441 (50.8648%)

Failed QC: 0 (0%)

Duplicates: 0 (0%)

Paired-end reads: 0 (0%)

100%マッピングされている。

samを評価する。

python RNAseqEval.py eval-mapping example_dataset/dmelanogaster_chr4_genome.fa mapped.sam -a example_dataset/dmelanogaster_chr4.gtf -o output
cat output #出力確認

user$ cat output

Command Line: RNAseqEval.py eval-mapping example_dataset/dmelanogaster_chr4_genome.fa mapped.sam -a example_dataset/dmelanogaster_chr4.gtf -o prefix

Reference format: ANNOTATION

General information:

Reference length = 1348131 bp

Number of chromosomes = 1

Chromosomes:

chr4: 5537446bp

Number of alignments in SAM file (total / unique) = 2833 / 2833

Alignments with / without CIGAR string = 2833 / 0

Mapping quality without zeroes (avg / min / max) = 44.56 / 1 / 47

Alignments with mapping quality (>0 / =0) = 5548 / 0