ロングリードを使った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
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
Number of matches / mismatches / inserts / deletes = 1186872 / 149935 / 34 / 7
Percentage of matches / mismatches / inserts / deletes = 0.89 / 0.11 / 0.00 / 0.00
Annotation statistics:
Total gene length = 4189315
Number of Transcripts / Exons (Multiexon transcripts) = 334 / 3131 (316)
Maximum number of exons in a gene = 46
Gene size (Min / Max / Avg) = 69 / 52075 / 12542.86
Exon size (Min / Max / Avg) = 6 / 9812 / 426.99
Mapping quality information:
Bases in reads (aligned / total) (percent) = (1336841 / 1336841) (100.00%)
Transcripts covered / missed / total = 230 / 104 / 334
Exons covered / missed / total = 1198 / 1933 / 3131
Alignments on transcript hit / missed = 2826 / 7
Alignments on exons hit / missed = 2826 / 7
Alignments hitting an exon (start / end / both) = 2939 / 2992 / 1590
Contiguous / non contiguous alignments: 1319 (46.54%) / 1514 (53.42%)
chr4(5537446 bp)へのアライメント結果のstatisticsがまとめられている。
詳細はGithubの解説を確認。
上ではbbmapを使いましたが、論文中ではBLASRやSTARのロングリード対応版などが使われています。興味がある人は論文を確認してください。
引用
Evaluation of tools for long read RNA-seq splice-aware alignment
Križanovic K, Echchiki A, Roux J, Šikic M
Bioinformatics. 2017 Oct 23
ロングリードのアライメント