DESeq2에서 연속적인 값을 condtion으로 받기

DESeq2에서는 2가지 이상의 condition에 대한 차이를 보고 Differentially Expressed Genes를 구할 수도 있지만 시간이나 약물 투여 농도에 따른 DEGs도 계산할 수 있다.

아래의 예시처럼 약물 농도는 0nM, 0.1nM, 1nM을 투여하고 이에 따른 유전자의 변화를 가지고 p-value를 구하게 된다.

Reference -

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

MultiQC 설치 및 실행하기

로그파일의 이름과 무관하게 형식만 맞으면 다 읽는 것으로 보이며 그럼에도 불구하고 실행시간은 오래걸리지 않는다. 

로그파일들의 크기가 대체로 작긴 하지만 폴더내에 있는 파일 중에 매치가되지 않으면 바로바로 점프하는 듯 하다.

설치 방법은 가장 간단한 방법은 PyPI을 사용하는 방법이다.

프로그램이 Python 기반이기 때문에 Python을 설치하고 PyPI를 사용하도록 하자.

2017/08/16 - [programming language/python] - Python 설치 및 실행하기

pip가 설치되어 있다면 아래처럼 입력하고

권한문제가 생긴다면 --user 옵션을 넣도록 하자.

--user 옵션을 사용했다면 실행 파일은 ~/.local/bin/multiqc 에 생겼을 것이다.

설치가 정상적으로 끝났다면 프로젝트 폴더에가서

명령을 수행하면 multiqc_report.html 파일이 생성된다.

html파일을 열면 프로젝트 폴더 안에 있는 대부분의 로그파일들이 올라와있는 것을 확인할 수 있다.

당연하게도 모든 중간작업마다 로그파일을 생성했어야만 정상적으로 파일을 읽고 작동한다.

table이나 plot 복사를 하면 report를 작성하는데 매우 유용하게 사용할 수 있을 것이다.

Reference-

http://multiqc.info/docs/#using-multiqc

RPKM, FPKM and TPM의 정의

RPKM (Reads Per Kilobase of transcript per Million mapped reads)

FPKM (Fragments Per Kilobase of transcripts per Million mapped reads)

TPM (Transcripts Per Million)

Reference -

http://www.incodom.kr/Expression_profiling

https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/

Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory in biosciences. 2012 Dec 1;131(4):281-5.

DESeq2 에서 multiple condition 수행하기

쭉 진행하면 나오는 결과는 test와 brain간의 비교이다. DESeq2는 default로 가장 처음과 마지막에 있는 컨디션을 비교하기 때문이다.

이를 아래처럼 수정하면 test와 UHR간의 비교 결과를 보여준다.

2, 3번째 값을 수정하면 그에 맞게 바뀌는데 당연히 이미 설정해놓은 컨디션 값을 지정해주어야하며

순서를 바꾸면 log2FoldChange값이 +에서 -로 바뀌게 된다.

log2FC값 외에 변경되는 값은 없다.

리눅스에서 프록시 설정하기

Centos 6에서 R 설치를 위한 라이브러리 설치

zlib

bzip

liblzma

pcre

curl

모든 dependency가 설치되었으면 R의 configure를 해보도록 하자.

중간에 warning이 뜨긴 했지만 일단 무시하고 진행했으며 R일 실행되는 것 까지 확인하였다.

Reference -

https://unix.stackexchange.com/questions/343452/how-to-install-r-3-3-1-in-my-own-directory

Arguments in R

argument를 더 복잡하게 쓰려면 optparse라는 라이브러리를 써도 되지만 간단하게 정리하고 싶다면 위와 같이 작성할 수 있다.

추가로 argument를 입력하지 않았을 때 간단한 설명을 넣고 싶다면 아래와 같이 하면 된다.

argument가 하나도 들어오지 않았다면 ERROR 메세지 뒤에 정해놓은 문자열을 출력하고 자동 종료된다.

Remove duplicates

source -