macでインフォマティクス

macでインフォマティクス

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

RNA seqのクオリティコントロールツール RSeQC

 

RSeQCはクオリティ、GCバイアス、PCRバイアス、ヌクレオチド組成バイアス、シーケンスのデプス、strandの特異性、カバレッジ均一性およびゲノムのfeature上のカバレッジ分布など、RNA-seq実験を総合的に評価するパッケージ。SAMとBAMを入力として利用できる。可視化にはRを利用している。

 

マニュアル

http://rseqc.sourceforge.net

 

インストール

依存

  • gcc
  • python2.7
  • numpy
  • R

pipで本体は導入できる。

pip install RSeQC

 

ラン

マニュアルからテストのbamやbed,gtfなどをダウンロードできる。ここではこのデータを使って動作を確認する。

 

ペアードエンドのbamからfastqに変換。

bam2fq.py -i test_PairedEnd_StrandSpecific_hg19.sam -o bam2fq_out1
  • -i    INPUT_FILE
  • -o   OUTPUT_PREFIX
  • -s   single-end for single-end sequencing. 

シングルエンドは-sをつける。

 

bamからwiggleフォーマット(UCSC解説)に変換。

UCSCのwigToBigWigコマンドが見つかれば、出力は自動でbigwigフォーマット(UCSC解説)に変換される。オーサーが準備したscript (fetchChromSizes) をマニュアルに従い準備しておく。

#bamはsamtools sortでcoordinate sortされている必要がある。
samtools sort -m 1000000000 input.bam input.sorted.bam
samtools index input.sorted.bam

#fetchChromSizesでスクリプトでサイズを計算
fetchChromSizes hg19 >hg19.chrom.sizes
bam2wig.py -s hg19.chrom.sizes -i sample.bam -o out -u
  • -s    Chromosome size file. Tab or space separated text file with 2 columns: first column is chromosome name/ID, second column is chromosome size. 
  • -u    Skip non-unique hit reads.
  • -t    Specified wigsum. Eg: 1,000,000,000 equals to coverage of 10 million 100nt reads. Ignore this option to disable normalization
  • -q    Minimum mapping quality for an alignment to be called “uniquely mapped”. default=30
  • -d   How read(s) were stranded during sequencing. For example: –strand=‘1++,1–,2+-,2-+’ means that this is a pair-end, strand-specific RNA-seq data, and the strand rule is: read1 mapped to ‘+’ => parental gene on ‘+’; read1 mapped to ‘-‘ => parental gene on ‘-‘; read2 mapped to ‘+’ => parental gene on ‘-‘; read2 mapped to ‘-‘ => parental gene on ‘+’. If you are not sure about the strand rule, run ‘infer_experiment.py’ default=none (Not a strand specific RNA-seq data).

 

bam/samのstatistics出力。

bam_stat.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam 
  • -q    Minimum mapping quality to determine uniquely mapped read.

出力。

$ bam_stat.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam

Load BAM file ...  Done

 

#==================================================

#All numbers are READ count

#==================================================

 

Total records:                          44826454

 

QC failed:                              0

Optical/PCR duplicate:                  0

Non primary hits                        0

Unmapped reads:                         0

mapq < mapq_cut (non-unique):           0

 

mapq >= mapq_cut (unique):              44826454

Read-1:                                 28674562

Read-2:                                 16151892

Reads map to '+':                       23884437

Reads map to '-':                       20942017

Non-splice reads:                       39959172

Splice reads:                           4867282

Reads mapped in proper pairs:           31803590

Proper-paired reads map to different chrom:0

 

 

リードのクリップアライメントの分布。

(clipped mappingが可能なアライナーを使っている必要がある)。

clipping_profile.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam -s "PE" -o out
  • -s    Sequencing layout. “SE”(single-end) or “PE”(pair-end).
  • -q    Minimum mapping quality (phred scaled) for an alignment to be considered as “uniquely mapped”. default=30

リード1の出力。

f:id:kazumaxneo:20180227225932j:plain

リード2も出力される。

 

 リードのmismatch分布。

mismatch_profile.py -l 101 -i sample.bam -o out
  • -l    Alignment length of read. It is usually set to the orignial read length. For example, all these cigar strings (“101M”, “68M140N33M”, “53M1D48M”) suggest the read alignment length is 101. [required]

 

リードのdeletion分布。

deletion_profile.py -i sample.bam -l 101 -o out
  • -l     Alignment length of read. It is usually set to the orignial read length.

 

 リードのinsertion分布。

insertion_profile.py -s "PE" -i test.bam -o out 
  • -s    Sequencing layout. “SE”(single-end) or “PE”(pair-end).

f:id:kazumaxneo:20180227234729j:plain

 

ランダムにリードをサンプリングしてbamをn個に分割

(ランダムだが均一なリード数になる)。

divide_bam.py -n 3 -i SingleEnd_StrandSpecific_50mer_Human_hg19.bam -o output
  • -n    Number of small BAM files

 

遺伝子上のカバレッジ分布

geneBody_coverage.py -r hg19.housekeeping.bed -i test.bam -o output
  • -r    Reference gene model in bed format. [required]

複数bamなら-i test1.bam,test2.bam,test3.bam 、またはbamのディレクトリを-i /data/ という風に指定する。bigwigから計算するgeneBody_coverage2.pyもある。

f:id:kazumaxneo:20180227233520j:plain

 

 

ペアードエンドの内側のインサートサイズを推測。

f:id:kazumaxneo:20180227234231j:plain

inner_distance.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam -o output -r hg19.refseq.bed12

f:id:kazumaxneo:20180301124426j:plain

 

スプライシングジャンクションをリファレンスと比較し、比率を返す。

junction_annotation.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam -o output -r hg19.refseq.bed12
  • f:id:kazumaxneo:20180228200527j:plain

f:id:kazumaxneo:20180227235016j:plain

 

 

レアなalternative splicingを探すためにデプスが十分か推測。

junction_saturation.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam -r hg19.refseq.bed12 -o output
  • -r    Reference gene model in BED format.

f:id:kazumaxneo:20180227235353j:plain

 

feature上のリードの分布

read_distribution.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam -r hg19.refseq.bed12

出力の定義はマニュアルで確認してください。

f:id:kazumaxneo:20180301124530j:plain

 

duplicateの分布。

read_duplication.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam -o output

 

 GC含量の分布

read_GC.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam -o output

  

f:id:kazumaxneo:20180228202855j:plain

リードの各部位の塩基の分布

read_NVC.py -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam -o output

 

クオリティの分布

read_quality.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam -o output_quality

f:id:kazumaxneo:20180301001942j:plain

f:id:kazumaxneo:20180301001944j:plain

 

各transcriptsごとのリードサイズ

RNA_fragment_size.py -r hg19.RefSeq.union.bed -i SRR873822_RIN10.bam > SRR873822_RIN10.fragSize

 

bamの分割(rRNA除去など)

split_bam.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam -r hg19.rRNA.bed -o output

 以下の3種類に分割される。rRNAをbedで与えれば、rRNAにマッピングされなかったリードをbamのまま抽出できる。

  1. XXX.in.bam: reads that are mapped to exon regions of the gene list (or reads consumed by gene list).
  2. XXX.ex.bam: reads that cannot be mapped the exon regions of the original gene list.
  3. XXX.junk.bam: qcfailed reads or unmapped reads.

 

ペアエンドbamからシングルエンドのbam2つに分割。

split_paired_bam.py -i Pairend_StrandSpecific_51mer_Human_hg19.bam -o output

 

TIN値からサンプルRNAの完全さを評価。

tin.py -r hg19.RefSeq.union.bed -i Pairend_StrandSpecific_51mer_Human_hg19.bam -o output

RIN値(1≦RIN≦10)は18Sおよび28SリボソームRNAの量に大きく依存したリボソームRNA完全性の尺度で、mRNA品質の直接測定ではない。実験ではリボゼロなどでrRNAを除いているため、下流のフローでRNAの完全さを分析できない。実際の状況は、「AU-リッチ配列」、「転写長」、「GC含量」、「二次構造」および「RNA-タンパク質複合体」などの因子に応じて、分解速度は転写物間で著しく異なることがある(ref:http://www.illumina.com/documents/products/technotes/technote-truseq-rna-access.pdf)。 これらの限界を克服するために、転写レベルでRNAの完全性を測定できるアルゴリズムTIN(リンク)が考案された。 TINは、発現された各転写物のスコア(0≦TIN≦100)を計算する。また全transcriptsの平均のmedTINはサンプル間のRNA完全性を比較するために使用できる。

 

 

ほかにも定量してFPKMを出すFPKM_count.pyや、RNA seqの種類を推測するinfer_experiment.py (strand specificかなどをマッピングから推測)、bigwigをノーマライズするnormalize_bigwig.py、2つのbigwigファイルを操作するoverlay_bigwig.py、ヘキサマーの頻度をリファレンスとリードで調べるread_hexamer.py、RPKMのがある。

 

引用

RSeQC: quality control of RNA-seq experiments.

Wang L, Wang S, Li W.

Bioinformatics. 2012 Aug 15;28(16):2184-5.