Illumina 시퀀싱에서 약 10-12개의 염기가 균등하게 분포하지 않는 패턴을 보인다. gDNA에서는 조금 더 드물지만 mRNA-seq에서는 대부분의 데이터가 이러한 패턴을 보이는데 이유를 찾아보았다.

Illumina에서는 이러한 현상의 원인을 랜덤프라이머를 제작하여 시퀀싱을 진행하지만 랜덤 프라이머가 완전한 랜덤이 아니기때문이라고 얘기한다.  

아래 plot은 bisulfite-seq 이다. bisulfite 처리로 인해 CtoT 변화로 C의 비율은 낮고 T의 비율이 높게 나온다. 하지만 그와 별도로 여전히 10개의 염기의 비율이 특이적이다.

해당 부분은 분석에 크게 영향을 주지 않으니 무시하고 진행하여도 상관없다.

Reference -

http://seqanswers.com/forums/showthread.php?t=11843

http://nar.oxfordjournals.org/content/38/12/e131

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/

RNA-seq 라이브러리 종류와 구별법

Susan M. Corely et al., Differentially expressed genes from RNA-Seq and functional enrichment results are affected by the choice of single-end versus paired-end reads and stranded versus non-stranded protocols, BMC Genomics, 2016

Single-end 

- single-end는 우리가 생각할 수 있는 가장 기본적인 형태이다. 단순히 길이가 약 35~150bp 정도가 되는 단일 reads들의 집합으로 존재하며 수 년 전에 생산한 데이터들은 이 형태를 많이 띄고있다. 다른 방법들에 비해 만들기 쉽지만 상대적으로 정확한 위치에 맵핑되지 않는다는 단점이 있기 때문에 Illumina사의 HiSeq등의 기술로 생산하는 short reads에서는 많이 보이지 않는 추세이다.

- 하지만 PacBio, Minion, MiSeq등의 기술은 reads의 길이를 350~수 kb에 해당하게 늘릴 수 있고 길이가 긴 reads는 정확환 위치에 맵핑할 확률이 늘어나기 때문에 특정한 상황에서는 여전히 쓰이고 있다.

Paired-end

- paired-end는 초기 single-end의 단점을 극복하고자 read를 양 쪽에서 읽어 한 쌍의 fragement를 만들어 더 정확한 위치에 맵핑되고자 하는 방식이다. 

Zhernakova, et al., PLoS Genet. 2013

- 위의 사진을 비교해보면 길이가 100bp reads가 있다고 가정할 때, single-end에서는 전체 genome상에 reads와 유사한 서열이 100bp만 있다면 그 곳에 맵핑 될 수 있지만 paired-end같은 경우에는 한 쌍이라서 100bp+약간의 간격+100bp가 모두 일치해야만 되기 때문에 200bp+@가 매치되어야한다. 이는 확률적으로 잘못된 위치에 맵핑될 확률이 single-end보다 낮아진다는 것을 의미한다. (@는 insert size를 말하며 라이브러리 제작시 원하는 만큼 줄 수 있다.)

Single-end와 Paired-end 구별법

- single-end와 paired-end의 구분법은 대부분의 paired-end는 왼쪽과 오른쪽 파일이 각각 존재한다. sample_1.fq와 sample_2.fq가 있다면 각 파일 안에 있는 readID가 같은 것 끼리 짝을 이루는 식이다. 한 파일 양 reads가 동시에 존재하는 경우도 있는데 이때도 파일 안에 readID를 자세히 보면 readID가 동일하고 끝에 1/2 이런식으로 구분되어 있다면 paired-end인 것이다.

Unstranded

- 라이브러리를 제작할 때 시퀀싱 기계가 reads를 읽기 위해서 양 말단에 adapter sequence를 붙이게 되는데 같은 서열의 adapter sequence가 붙게 되면 reads의 앞 뒤를 구분할 수 없게 된다. 

Strand-specific

- Strand-specific은 양 말단의 adapter sequence를 다르게 붙임으로써 reads의 앞 뒤를 구분해 주는 방식을 말하는데, 제작방식은 크게 세 가지가 존재하며, 각 방식에따라 맵핑할 때 옵션을 다르게 줘야 할 수도 있으니 구분할 수 있는것이 좋다.

- 맵핑 프로그램을 보면 fr-strand, rf-strand 등이 존재하는데 앞의 fr은 foward reverse, rf는 reverse foward를 의미하며 각각 앞 뒤의 reads가 정방향인지 역방향인지를 의미한다. 어떤 것을 선택할 지는 paired-end 제작 방식에 따라 달라질 수 있기 때문에 지정해주어야한다. (STAR는 프로그램이 알아서 계산하지만 tophat은 직접 넣어줘야 한다. 아래는 tophat의 메뉴얼중의 일부로 옵션별 설명이 있어서 캡쳐해서 가져와 보았다.)

- 아래는 strand 예측이 왜 중요한지를 보여주는 figure이다. 가끔 genome 상에서 유전자들이 겹치는 부분이 존재한다. (DNA상에서의 방향성이 같고 겹친다면 같은 유전자이기 때문에 당연히도 예시의 두 유전자의 방향성은 반대이다.) Unstranded RNA-seq signal을 보면 어떤 유전자에서 나온 reads인지 구분이 안되지만 stranded RNA-seq이라면 RPD signal처럼 구분할 수 있다. 또한 발현량 분석 프로그램이 유전자의 방향성과 반대의 방향성을 가진 reads들은 유전자의 발현량으로 계산하지 않는 다는 점에서 strand-specific은 더 정확한 유전자 발현량을 계산할 수 있게 해주는 방식이라고 말할 수 있다.

(해당 논문은 unstrand데이터를 stranded데이터로 바꾸는 프로그램을 소개하며 해당 방식을 적용했을 때 민감도와 특이도가 높아지는 것을 보여준다.)

Bo-Hyun You et al., High-confidence coding and noncoding trnascriptome maps, Genome Research, 2017

Unstranded와 Strand-specific 구별법

- 아래 그림을 보면 방향성이 <<< 인 유전자에 있는 reads는 대부분 F1R2를 가지고 있고 >>>에는 F2R1를 가지고 있는 것을 확인할 수 있다. (가끔 by chance로 아닌 reads들이 있지만 5% 미만이다.)

위의 방식을 지원하는 프로그램도 존재한다.

2018/08/06 - [bioinformatics] - RSeQC를 사용하여 stranded 데이터 확인하기

마지막으로 Single-end, Paired-end와 Stranded, Unstranded는 별개의 개념이라는 것을 알아야 한다. Single-end이더라도 Stranded일 수도 있고 Unstranded일 수도 있으며 Paired-end도 마찬가지이다.

Reference -

https://genome.cshlp.org/content/27/6/1050.full?sid=63335fd5-50fe-46bd-b849-c88bbc2a0a84

https://www.biostars.org/p/70833/#70834

http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1003594