RSeQC를 사용하여 stranded 데이터 확인하기

지난 글에 이어 RNA-seq에서 데이터가 stranded인지 unstranded인지 구별하는 방법에 대해 더 알아보기로 한다.

2018/08/03 - [bioinformatics] - RNA-seq 라이브러리 종류와 구별법

지난 글에서는 IGV로 직접 reads를 확인하는 방법을 설명했지만 이를 확인시켜주는 프로그램도 당연히 존재한다.

하지만 개인적으로 RNA-seq데이터를 다룰 때는 한 번쯤은 IGV로 데이터가 분석에 사용하기에 적절한 품질을 지녔는지 확인해보는것이 좋다고 생각하기 때문에 직접 확인하는 것을 권장하지만 데이터가 수십, 수백 개가 된다면 일일히 다 확인하는 것은 불가능하다.

그럴때 RSeQC에 포함되어있는 infer_experiment.py 스크립트를 사용하여 stranded인지 여부를 확인해보도록 하자.

#http://rseqc.sourceforge.net/#infer-experiment-py

위의 링크를 타고 들어가면 자세한 설명이 있기때문에 여기선 간단하게 요약만 한다.

필요한 데이터는 Alignement file (BAM,SAM) 파일과 annotation.bed 파일이다.

전자는 당연히 준비가 되어 있을 것이고 후자는 RSeQC 홈페이지에서 아래 링크로 가면 받을 수 있다. (human과 mouse만을 지원한다. 다른 종은 직접 만들거나 다른 데이터베이스에서 찾기 어렵지 않을 것이다.)

http://rseqc.sourceforge.net/#download-gene-models-update-on-08-07-2014

결과해석은 위의 infer-experiment-py 링크에 설명을 참고하면 된다.

infer_experiment.py -r hg19.refseq.bed12 -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam

#Output::

This is PairEnd Data
Fraction of reads failed to determine: 0.0172
Fraction of reads explained by "1++,1--,2+-,2-+": 0.4903
Fraction of reads explained by "1+-,1-+,2++,2--": 0.4925

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

#Output::

This is PairEnd Data
Fraction of reads failed to determine: 0.0072
Fraction of reads explained by "1++,1--,2+-,2-+": 0.9441
Fraction of reads explained by "1+-,1-+,2++,2--": 0.0487

infer_experiment.py -r hg19.refseq.bed12 -i SingleEnd_StrandSpecific_36mer_Human_hg19.bam

#Output::

This is SingleEnd Data
Fraction of reads failed to determine: 0.0170
Fraction of reads explained by "++,--": 0.9669
Fraction of reads explained by "+-,-+": 0.0161

단순하게 위의 Fraction과 아래의 Fraction이 약 50:50을 이룬다면 annotation되어 있는 유전자의 방향성과 맵핑된 reads의 방향성간의 통일성이 없기 때문에 unstranded 데이터라고 말할 수 있으며 반대로 첫 번째 Fraction과 두 번째 Fraction중에 한 곳으로 몰려있다면 방향성이 통일되어 있다는 뜻이기 때문에 stranded라고 할 수 있다.

위의 Fraction과 아래 Fraction중에 어디에 몰리는지는 stranded 데이터를 만들 때 방법에 차이에서부터 기인한다. library 설명 포스팅을 참고하면 된다.

Reference -

http://rseqc.sourceforge.net/