macでインフォマティクス

macでインフォマティクス

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

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

本体 Github

https://github.com/kkrizanovic/RNAseqEval

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 K1, Echchiki A2,3, Roux J2,3,4, Šikic M1,5. 

Bioinformatics. 2017 Oct 23. doi: 10.1093/bioinformatics/btx668.

 

 

ロングリードのアライメント

https://github.com/PacificBiosciences/cDNA_primer/wiki/Aligner-tutorial%3A-GMAP%2C-STAR%2C-BLAT%2C-and-BLASR