반응형

RPKM, FPKM and TPM의 정의




RPKM, FPKM, TPM은 생물정보학에서 상당히 쉽게 접할 수 있는 용어들이다.

여러 샘플의 RNA-seq을 발현 분석할 때 정규화(normalize)된 발현량을 의미하는 것으로 과거에는 RPKM이 많이 쓰였으나 점차 FPKM을 거쳐 TPM을 많이 사용하고 있다.

각 용어가 의미하는 바는 크게 다르지 않으며 샘플간의 차이 (주로 total depth)로 인한 발현량을 정규화하여 비교하기 위해 사용하게 되었다.

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

- 추후에 FPKM 설명에서 FPKM과 RPKM에 차이에 대해서 설명하겠지만 과거에 single-end RNA-seq이 주로 생산되었을 때 단순히 read counts in transcripts, gene length, total mapped reads 만을 가지고 계산한 값이다.

enter image description here

- 먼저 total mapped reads로 나누는 이유는 sample에 따라서 total depth가 다른 것을 normalize해주기 위함이다. 
sample A는 1,000,000개의 reads가 있고 sample B에는 2,000,000개의 reads가 있을 때, 같은 유전자에 같은 갯수의 reads가 붙었다고하면 sample A가 더 발현량이 높다고 계산하기 위함이다. 
(sample B가 PCR cycle이 한 번 더 진행되었기 때문에 전체 depth가 높은 것이라고 가정한 것이다. 단순히 reads count가 같다고 해서 발현량이 같은 것이 아니다!)

- gene length로 나누는 이유는 샘플간의 비교가 아닌 유전자 간의 비교를 위함이다. 
gene A는 길이가 1000고 gene B는 길이가 2000일 때 위와 마찬가지로 같은 갯수의 reads가 붙었다면 gene A의 발현량이 더 높다고 추측할 수 있다. 
(길이가 2배이고 발현량이 같다면 gene B가 두 배 많은 reads를 생산 되었어야 한다!)

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

- FPKM과 RPKM의 차이는 paired-end로 생산된 RNA-seq에서 나타난다. paired-end는 하나의 reads에서 두 개의 fragments가 나온다고 생각하면 된다. 즉 left fragements와 right fragments가 각각 계산되어 RPKM의 대략 두 배의 값을 가지게 된다.(paired-end의 두 fragments가 항상 같이 맵핑되는 것은 아니기 때문에 정확하게 두 배 일 수는 없다.) 
용어가 헷갈릴 수 있는데 일반적으로 paired-end는 reads를 두 개 가지고 있다고 말하는 경우도 있다. 하지만 여기서는 각각을 fragments로 계산하고 있음을 유의해야 한다.


TPM (Transcripts Per Million)

- TPM에서는 total mapped reads를 사용하지 않으며 transcripts level에서의 값이 계산된다. gene level로 계산하고 싶으면 TPM값들을 모두 더하면 된다.

계산하는 방법은 순서가 중요한데

1. 각각의 transcipt에 대해 mapped reads / transcripts length / killobase 로 값을 구한다. (=normalized transcripts expression)
2. 1번에서 구한 값들을 모두 더한뒤 1,000,000으로 나눈다. (=scaling factor)
3. 각각 1번에서 구한 값들을 2번에서 구한 값으로 나눈다. (=TPM)


아래 reference에 있는 논문에 따르면 RPKM을 사용하는 것 보다 TPM을 사용하는 것이 더 정확한 발현 값을 구할 수 있다고 한다.


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.


반응형

'bioinformatics' 카테고리의 다른 글

RNA-seq 라이브러리 종류와 구별법  (1) 2018.08.03
MultiQC 설치 및 실행하기  (1) 2018.07.31
Remove duplicates  (0) 2018.07.24
Gene ontology analysis - DAVID  (0) 2018.07.16
_PAR_Y in Genecode annotation  (0) 2018.07.13
반응형

DESeq2 에서 multiple condition 수행하기




DESeq2는 여러 컨디션이 있어도 DEG분석이란 결국 pairwise하게 비교해야 하기 때문에 condtion을 구분 해 놓은뒤 결과를 원하는 두 컨디션을 놓고 비교해서 각각 분석해야 한다.

SRA010153 데이터로 전체적인 DEG analysis workflow는 추후에 정리하기로 하고 여기서는 multiple condition의 비교에만 중점을 두고 설명하겠다.

SRA010153은 여러 컨디션이 있지만 Brain과 UHR phi X 데이터만 받은 상태에서 임의로 test 컨디션을 아래처럼 추가했다. (원래는 7개의 Brain과 7개의 UHR이다.)

sampleCondition<-c('Brain','Brain','Brain','Brain','Brain','Brain','test', 'UHR','UHR','UHR','UHR','test','UHR','test')

sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition)    

                                                                   

ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~condition)

write.csv(assay(ddsHTSeq),quote=F,file="raw_count.csv")                                                                               

                                                

colData(ddsHTSeq)$condition<-factor(colData(ddsHTSeq)$condition, levels=c('Brain','UHR','test'))

                                                                                        

cts <- counts(ddsHTSeq)                                                                 

geoMeans <- apply(cts, 1, function(row) if (all(row == 0)) 0 else exp(mean(log(row[row != 0]))))

ddsHTSeq <- estimateSizeFactors(ddsHTSeq, geoMeans=geoMeans)                            

dds<-DESeq(ddsHTSeq)


res <- results(dds)


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


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


results(dds,contrast=c('condition','test','UHR'))


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


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


> results(dds,contrast=c('condition','UHR','test'))

log2 fold change (MLE): condition UHR vs test 

Wald test p-value: condition UHR vs test 

DataFrame with 60461 rows and 6 columns

                          baseMean     log2FoldChange             lfcSE

                         <numeric>          <numeric>         <numeric>

ENSG00000000003   92.1871316897525  0.479423175690374  0.43657650183317

ENSG00000000005   2.33464385903914  0.243364522893621 0.807778836598724


> results(dds,contrast=c('condition','test','UHR'))

log2 fold change (MLE): condition test vs UHR 

Wald test p-value: condition test vs UHR 

DataFrame with 60461 rows and 6 columns

                          baseMean      log2FoldChange             lfcSE

                         <numeric>           <numeric>         <numeric>

ENSG00000000003   92.1871316897525  -0.479423175690374  0.43657650183317

ENSG00000000005   2.33464385903914  -0.243364522893621 0.807778836598724


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





반응형
반응형

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




외부로 연결되는 서버 A가 있고 이 서버를 통해서 연결할 수 있는 내부만 연결이 되어 있는 서버 B에서 작업할 때 프록시 설정을 해 놓으면 A서버가 중개해서 B서버가 외부와 통신할 수 있게 된다.

B서버의 .bash_profile 안에 아래와 같이 입력하면 된다.

#proxy setting
export http_proxy=http://[ID]:[PASSWORD]@ip:port/
export https_proxy=http://[ID]:[PASSWORD]@ip:port/
export ftp_proxy=http://[ID]:[PASSWORD]@ip:port/

ID와 password가 A서버와 B서버가 동일하다면 생략해도 되며 ip와 port는 반드시 넣어주어야 한다.

port는 proxy의 기본 port number가 3128라서 따로 설정하지 않았다면 3128을 넣으면 된다.



반응형
반응형

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




Root권한이 있다면 yum이나 apt-get으로 설치하면서 진행할 수 있겠지만 local로 설치하려고 하니 local library에 따로 설치해야 하는 것들이 많아서 정리해 보고자 한다.

먼저 local library 설치 폴더를 만든다. 간단하게 만들어도 무방하다.

mkdir ~/library 

zlib

제일 처음 문제가 생긴 부분은 zlib 버전이 1.2.5 미만이라는 것. zlib를 다운받는다. 이때 1.2.10 이상 버전을 받으면 lexiographically 1.2.10이 더 낮다고 판단하는 문제가 있기 때문에 무난하게 1.2.9로 받도록 하자. (lexiographically는 숫자를 문자로 인식해서 5와 15를 비교할 때 5 > 1 로 비교해서 5가 15보다 더 높다고 비교하는 방식이다.)
참조 (https://stackoverflow.com/questions/42076936/zlib-bz2-library-and-headers-are-requried-for-compiling-r)

wget https://osdn.net/frs/g_redir.php?m=kent&f=libpng%2Fzlib%2F1.2.9%2Fzlib-1.2.9.tar.gz
tar -zxf zlib-1.2.9.tar.gz
cd zlib-1.2.9
./configure --prefix=~/library
make && make install

일반적인 설치 방법대로 하면 zlib가 설치된다. 설치가 끝나면 ~/library 폴더 안에 lib와 include폴더가 생성되는데 여기에 나머지 라이브러리들도 다 넣고 이 라이브러리 경로를 잡아서 R을 설치 할 것이다.

bzip

bzip은 설치방법이 약간 다르다. configure가 없다. 바로 make를 하는데 prefix가 대문자임에 주의하면서 진행하면 된다.

wget http://www.bzip.org/1.0.6/bzip2-1.0.6.tar.gz
tar -zxf bzip2-1.0.6.tar.gz
cd bzip2-1.0.6
make --PREFIX=~/library

liblzma

lzma는 lzma를 직접 설치하는게 아니라 XZ를 설치하면 자동으로 해결된다. https://tukaani.org/xz/ 에서 받을 수 있다.

wget https://tukaani.org/xz/xz-5.2.4.tar.gz
tar -zxf xz-5.2.4.tar.gz
cd xz-5.2.4
./configure --prefix=~/library
make && make install

pcre

pcre는 추가로 설정해줘야하는 부분이 있다. 우선 파일은 https://ftp.pcre.org/pub/pcre/에서 받으면 되는데 pcre2를 받지 말고 pcre 8.42정도를 받자.

wget https://ftp.pcre.org/pub/pcre/pcre-8.42.tar.gz
tar -zxf pcre-8.42.tar.gz
cd pcre-8.42
./configure --prefix=~/library --enable-utf8
make && make install

cd ~/library/include
mkdir pcre
ln -s pcre* include/.

make 이후에 ~/library/include에 생성된 pcre.h 헤더 파일이 있는데 이 파일을 include/pcre 폴더 안에 링크로 넣어주는 작업까지 해야 완료가 된다. 왜인지는 모르겠지만 R 설치시 헤더파일을 pcre폴더 안에서 찾기 때문이다.

curl

마지막으로 curl은 일반적인 설치 방식을 따라가면 된다. https://curl.haxx.se/download.html

wget https://curl.haxx.se/download/curl-7.61.0.tar.gz
tar -zxf curl-7.61.0.tar.gz
cd curl-7.6.10
./configure --prefix=~/library
make && make install



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


./configure --prefix=/PATH/TO/INSTALL/R --enable-R-shlib LDFLAGS="-L~/library/lib/" CPPFLAGS="-I~/library/include/"

make && make install




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


Reference -

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

반응형

'Computer Science > R' 카테고리의 다른 글

Kegg pathway에 속하는 유전자 정보 가져오기  (2) 2018.09.11
pheatmap으로 heatmap그리기  (0) 2018.09.11
R에서 Dataframe 합치기  (0) 2018.09.05
DESeq2 에서 multiple condition 수행하기  (1) 2018.07.27
Arguments in R  (0) 2018.07.25
반응형

Arguments in R




Rscript를 사용할 때 argument를 input으로 받는 방법.

args = commandArgs(trailingOnly=TRUE)
species <- args[1]
inputfile <= args[2]

Rsciprt test.R human hg19.fasta


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


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


if(length(args)==0 {

        stop("All argument must be supplied ex) human hg19.fasta",call.=FALSE))

}


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



반응형

+ Recent posts