반응형

Gene id conversion in R




R에서 gene id로부터 다른형식의 geneid 값을 가져오는 방법에 대해서 설명하고자 한다.

예시는 human의 ensembl geneid를 입력값으로 받아 hgnc_symbol로 바꾸는 것이지만 종을 다르게 하거나 hgnc_symbol이 아닌 다른 정보도 얼마든지 가져올 수 있다.

users guide 주소 :


library(biomaRt)

mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol"),values=df$Geneid,mart= mart)

위의 코드는 df$Geneid에 ensembl geneid가 있는 상태에서 매칭되는 hgnc_symbol을 가져와 g_list에 저장한 것이다. 

useDataset에서 "hsapiens"를 다른 종으로 바꿀 수 있으며

getBM에서 attributes를 hgnc_symbol이 아니라 다른정보 (enterzgene, refseq_mrna, interpro, interpro_description 등)으로 바꾸면 해당 정보를 가져올 수 있다. 


물론 아래처럼 사용하여 동시에 가져올 수도 있다.


ipro = getBM(attributes=c("refseq_mrna","interpro","interpro_description"), filters="refseq_mrna", values=refseqids, mart=mart)


## refseq_mrna interpro interpro_description ## 1 NM_000546 IPR002117    p53 tumour suppressor family ## 2 NM_000546 IPR008967    p53-like transcription factor, DNA-binding ## 3 NM_000546 IPR010991    p53, tetramerisation domain ## 4 NM_000546 IPR011615    p53, DNA-binding domain ## 5 NM_000546 IPR012346    p53/RUNT-type transcription factor, DNA-binding domain superfamily


어떤 정보를 가져올 수 있는지는 위의 users guide를 참조하기 바란다.

listAttributes(mart)를 입력하면 사용한 데이터셋에 대한 가능한 attributes가 나온다.


> head(listAttributes(mart),10)

                            name                  description         page

1                ensembl_gene_id               Gene stable ID feature_page

2        ensembl_gene_id_version       Gene stable ID version feature_page

3          ensembl_transcript_id         Transcript stable ID feature_page

4  ensembl_transcript_id_version Transcript stable ID version feature_page

5             ensembl_peptide_id            Protein stable ID feature_page

6     ensembl_peptide_id_version    Protein stable ID version feature_page

7                ensembl_exon_id               Exon stable ID feature_page

8                    description             Gene description feature_page

9                chromosome_name     Chromosome/scaffold name feature_page

10                start_position              Gene start (bp) feature_page



G_list에 해당 정보를 담았다면 이를 기존의 df와 합치는 과정이 필요하다. 



1. df$Geneid와 G_list$ensembl_gene_id의 값이 같을 때 두 data frame을 합치는 방식이다.


원래는 이 방법을 사용하고 있었으나 ensembl_gene_id가 위의 데이터 베이스에 없을 때 결과 df의 사이즈가 입력할 때와 달라지는 것을 확인하여 2번의 방법을 사용하는 것을 추천한다.


df <- merge(df,G_list,by.x="Geneid",by.y="ensembl_gene_id")



2. df에 hgnc_symbol 열을 미리 만들고 내용은 공란으로 채워넣는다. 공란으로 채우는 이유는 ensembl gene id가 데이터 베이스 없거나 또는 ensembl gene id는 있지만 여기에 매칭되는 hgnc symbol이 없어도 행을 유지시키기 위함이다.


아래의 코드를 사용하면 df$Geneid와 G_list$ensembl_gene_id가 매칭될 때 G_list$hgnc_symbol의 값을 df$hgnc_symbol에  넣는 다는 의미이다.


df$hgnc_symbol = ""
df["hgnc_symbol"] = lapply("hgnc_symbol", function(x) G_list[[x]][match(df$Geneid, G_list$ensembl_gene_id)])

df$hgnc_symbol을 확인해보면 값이 없는 부분은 공란으로 남아있고 match된 부분은 모두 ensembl gene id에 대응하는hgnc_symbol값이 들어 있을 것이다.


Reference -

https://www.bioconductor.org/packages/devel/bioc/vignettes/biomaRt/inst/doc/biomaRt.html




반응형
저작자표시 비영리 변경금지

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

pheatmap 값에 따른 color 범위 조절하기  (0) 2018.09.12
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
반응형

pheatmap 값에 따른 color 범위 조절하기




heatmap에서는 일반적으로 발현량을 RdBl 색깔로 표현한다.


p-value를 여기서 적용하면 0.05라는 경계를 두고 보기에는 조금 불편하다. 0.6과 0.4를 사람 눈으로는 구별 할 수 없기 때문이다.


value에 따라서 구분하는 방법을 설명하고자 한다.




값이 -1부터 1까지 존재한다. ( DEG 분석에서 발현이 증가하는지 감소하는지 여부를 나타내려고 했다. )

1-adjust pvalue니까 0.95이상이거나 -0.95이하일 때가 p-value 0.05 이하이다.


library(RColorBrewer)


pheatmap(df, legend_breaks =  c(-0.95, 0, 0.95, max(df)), legend_labels = c("-0.95", "0", "0.95", "1-(p.adj)\n"), legend=T,  color = RColorBrewer::brewer.pal(7,"RdYlBu"), breaks = c(-1,-0.95,-0.9,-0.3,0.3,0.9,0.95,1))


legend_breaks = 범례에 글이 표시될 부분을 마크한다. 글자 크기가 너무 커서 전부를 표시하려 하지는 않았다.

legend_labels = breaks가 가리키는 포인트에 입력될 문자열을 넣으면 된다.

legend = T 범례를 표시할 것인지 여부. default는 F이다.

color = plot에 들어갈 색깔을 지정한다. 여기서는 Read/Yellow/Blue의 순서대로 7개로 나누어진 palette를 의미한다.

breaks = 데이터를 어떻게 나눌지를 의미한다. -1에서 -0.95 사이의 값이 color에서 1번 색깔로 지정 될 것이다.


breaks를 더 늘리게 되면 palette의 숫자 7도 여기에 맞춰서 늘려주면 된다. 색을 세 개를 혼합하였다면 가능하면 홀수로 맞추도록 하자.



Reference -

https://stackoverflow.com/questions/32545256/define-specific-value-colouring-with-pheatmap-in-r




반응형
저작자표시 비영리 변경금지

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

Gene id conversion in R  (0) 2018.11.16
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
반응형

Kegg pathway에 속하는 유전자 정보 가져오기



kegg pathway에는 하나에 pathway에 수십에서 수 백가지 유전자가 포함되어 있다. 

R에서 특정 pathway에 포함되어 있는 유전자 리스트를 불러와서 발현 값을 본다던지 하는 분석을 할 수 있는데

여기서 유전자 리스트를 불러오는 법을 설명하고자 한다.



(04010 pathway는 295개의 유전자를 포함하고 있다.)


아래 코드는 kegg pathway에 해당하는 ensembl gene id list를 가져오는 코드이다.


library(KEGGREST)


pathid = paste("path:hsa",i,sep="")

kegggenes = keggLink("hsa",pathid)

genelist = list()

for (j in 1:length(kegggenes)) {

ensemblgenes = keggGet(c(kegggenes[[j]]))[[1]]$DBLINKS

for (k in ensemblgenes) {

if (startsWith(k,"Ensembl:") == T) {

geneid = strsplit(k," ")[[1]][2]

genelist = append(genelist, geneid)

}

}

}


keggLink("hsa","path:hsa04010")를 입력하면 04010 pathway에 대한 모든 정보를 가져오게 되고 여기에 gene id는 DBLINKS에 있는데 이걸 다시 가져와야한다.

DBLINKS에는 ensembl 뿐만 아니라 여러 종류의 gene id를 모두 가지고 있기 때문에 startswith를 사용하여 Ensembl로 시작하는 gene id만을 genelist 라는 변수에 저장하였다.



반응형
저작자표시 비영리 변경금지

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

Gene id conversion in R  (0) 2018.11.16
pheatmap 값에 따른 color 범위 조절하기  (0) 2018.09.12
pheatmap으로 heatmap그리기  (0) 2018.09.11
R에서 Dataframe 합치기  (0) 2018.09.05
DESeq2 에서 multiple condition 수행하기  (1) 2018.07.27
반응형

pheatmap으로 heatmap그리기




pheatmap은 pretty heamaps의 약자로 heatmap을 그릴 때 더 이쁘고 쉽게 그려보자는 취지에서 만든 R 패키지이다. 주로 사용하게 되는 몇 가지 예시에 대해서 더 자세하게 정리해보고자 한다.


첫 스탭은 당연히 pheatmap을 불러오는 것.


library("pheatmap")



두 번째는 범례를 만드는 법이다. 

데이터 frame에 세포주와 약물을 처리했을 때 샘플들이 많으면 한 눈에 구분하기가 쉽지 않다. 구분하기 쉬운 색상을 미리 골라서 보기 쉽게 분리하도록 하자.

cell_line = c("cell_1","cell_2","cell_3","cell_4","cell_5")
drug = c("drug_1","drug_2","durg_3")

## annotation dataframe
anno.df = data.frame(cell_line=anno_cell_line, drug=anno_drug)
rownames(anno.df) = anno_samples

## annotation color      
ann_colors = list(
        cell_line = c("cell_1" = "#235725", cell_2 = "#1FD7F1", cell_3 = "#E4AF30", cell_4 = "#CADF7C", cell_5 = "#Fc0D7D"),
        drug = c(drug_1 = "#BDBE32", drug_2 = "#FE9089", "drug_3" = "#82B7FC")
)

위의 예시는 cell 5 종류와 drug 3 종류를 구분하였고 data frame을 만들었다. 
anno_samples는 cell_line과 drug 사이에 _를 넣고 붙인 리스트이다. 나중에 데이터와 일치화 시켜야 하기 때문에 상황에 맞게 다르게 줘야 할 수도 있다.

data frame이 제대로 만들어 졌다면 annotation 정보를 가지고 있는 아래와 같은 구조를 가질 것이다.

                cell_line       drug
cell_1_drug_1   cell_1  drug_1
cell_1_drug_2   cell_1  drug_2
cell_1_drug_3   cell_1  drug_3
cell_2_drug_1   cell_2  drug_1
cell_2_drug_2   cell_2  drug_2
cell_2_drug_3   cell_2  drug_3
cell_3_drug_1   cell_3  drug_1
cell_3_drug_2   cell_3  drug_2
cell_3_drug_3   cell_3  drug_3
cell_4_drug_1   cell_4  drug_1
cell_4_drug_2   cell_4  drug_2
cell_4_drug_3   cell_4  drug_3
cell_5_drug_1   cell_5  drug_1
cell_5_drug_2   cell_5  drug_2
cell_5_drug_3   cell_5  drug_3

이후에 데이터가 들어가 있는 data frame과 rowname을 일치시켜주고 plot을 그리면 된다.

names(df) = anno_samples


pdf("test.pdf")

pheatmap(df, cluster_rows=F, cluster_cols=F, show_rownames=T, annotation_col=anno.df, annotation_colors = ann_colors, main=("pheatmap_test"), fontsize_row=6, legend_breaks =  c(-0.5, 0, 0.5, max(newdf)), legend_labels = c("-0.5", "0", "0.5", "1-(q-value)\n"), legend=T)

dev.off()


＊여기서 annotation 정보가 있는 data frame의 row name이 데이터가 들어가 있는 data frame의 column name과 일치되어야 한다. 

rownames(anno.df) = anno_samples

names(df) = anno_samples


행과 열을 clustering하는 여부나 이름 표시 여부 등은 아래 메뉴얼을 참고하는 것이 더 빠를 듯 하다. 대부분은 직관적으로 이해할 수 있다.


legend_breaks와 legend_labels은 데이터 표시 범위의 범례를 조절하는데 default로 놓아도 크게 무리없는듯 하다.


예시대로 plot을 그리면 아래처럼 나온다. 색상은 직접 조절하도록 하자.


2018/08/31 - [etc.] - 16진수 RGB코드 알아내는법






Reference -

https://cran.r-project.org/web/packages/pheatmap/pheatmap.pdf



반응형
저작자표시 비영리 변경금지

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

pheatmap 값에 따른 color 범위 조절하기  (0) 2018.09.12
Kegg pathway에 속하는 유전자 정보 가져오기  (2) 2018.09.11
R에서 Dataframe 합치기  (0) 2018.09.05
DESeq2 에서 multiple condition 수행하기  (1) 2018.07.27
centos 6에서 R 설치를 위한 라이브러리 설치  (2) 2018.07.25
반응형

R에서 Dataframe 합치기




main이 같은 dataframe끼리 합친다.


df <- merge(df1, df2, by="main")


두 dataframe의 row name이 다르다면 각각을 지정해준다.


df <- merge(df1, df2, by.x="xmain", by.y="ymain")


df1에서는 xmain이라는 row와 df2에서는 ymain이라는 row의 값이 같으면 합친다.


값이 채워지지 않는다면 빈칸으로 존재하는데 이를 그냥 무시하고 지나가면 값이 밀릴 수 있다.


df$xmain <- ifelse(df$xmain == "" , "NA", df$xmain)


df의 xmain이 비어있다면 "NA"로 채우고 비어있지 않다면 값을 유지하고 지나간다.





반응형
저작자표시 비영리 변경금지

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

Kegg pathway에 속하는 유전자 정보 가져오기  (2) 2018.09.11
pheatmap으로 heatmap그리기  (0) 2018.09.11
DESeq2 에서 multiple condition 수행하기  (1) 2018.07.27
centos 6에서 R 설치를 위한 라이브러리 설치  (2) 2018.07.25
Arguments in R  (0) 2018.07.25
반응형

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값 외에 변경되는 값은 없다.





반응형
저작자표시 비영리 변경금지

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

Kegg pathway에 속하는 유전자 정보 가져오기  (2) 2018.09.11
pheatmap으로 heatmap그리기  (0) 2018.09.11
R에서 Dataframe 합치기  (0) 2018.09.05
centos 6에서 R 설치를 위한 라이브러리 설치  (2) 2018.07.25
Arguments in R  (0) 2018.07.25
반응형

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 메세지 뒤에 정해놓은 문자열을 출력하고 자동 종료된다.



반응형
저작자표시 비영리 변경금지

'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
centos 6에서 R 설치를 위한 라이브러리 설치  (2) 2018.07.25

+ Recent posts

티스토리툴바