반응형

python 파일 입출력

 

 

python에서 파일을 찾고, 파일이 없을시 except 출력, 있다면 파일을 열어서 본문을 수행하고 파일 닫고 종료.

 

        try : 
                list_file = open(list_file_name) 

        except OSError : 
                print('can\'t find '+list_file_name) 

        else : 
                with listfile : 
                        for line in listfile :
                        	# do somthing

 

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

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

Primer 서열 분석을 위한 python 코드  (0) 2021.08.17
String Format으로 길이 고정하기  (0) 2020.06.24
python multi-level argparse  (0) 2019.07.12
Python 설치 및 실행하기  (0) 2017.08.16
cannot mkdir R_TempDir 에러  (0) 2016.09.07
반응형

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
반응형

tar 디렉토리 지정해서 압축 풀기




원하는 경로에 압축을 풀고싶다면 아래처럼 입력하면 된다.


tar zxvf file.tar.gz -C path


해당 폴더에 파일을 바로 압축을 풀어주니 미리 빈 폴더를 생성하고 그 폴더 안에 푸는것이 좋을 것 같다.



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

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

비밀번호 없이 ssh로 로그인하기  (0) 2019.08.26
linux 계정 정보 옮기기  (1) 2019.08.21
특정 파일이나 디렉토리를 남기고 모두 지우기  (0) 2018.09.03
리눅스에서 프록시 설정하기  (0) 2018.07.26
neocomplcache vim plugin 설치하기  (0) 2018.07.12
반응형

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
반응형

Argument받아서 처리하기




bash 스크립트로 argument를 여러 개 받아 해당 수 만큼 for loop을 돌릴 때 쓰는 방법이다.


전체 argument를 반복하고 싶다면 아래처럼 하면 된다.


for i in $@ ; do echo $i ; done 


특정 구간만을 반복시키고 싶다면 indexing을 넣어주어야 한다. 아래 예시는 2번부터 끝까지 반복이다.


for i in ${@:2} ; do echo $i ; done


2번부터 3번 반복하고 싶다면 (즉 2부터 4까지) 아래처럼 더 추가한다.


for i in ${@:2:3} ; do echo $i ; done



argument가 몇 개 인지 알고싶다면 "$#" 를 사용하면 된다.


argument가 0개 일때의 조건문.


if [ $# -eq 0 ] ; then

        echo "Number of argument is 0"

else 

        echo "Number of argument is $#"

fi





Reference -

https://stackoverflow.com/questions/255898/how-to-iterate-over-arguments-in-a-bash-script

https://stackoverflow.com/questions/4423306/how-do-i-find-the-number-of-arguments-passed-to-a-bash-script

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

특정 파일이나 디렉토리를 남기고 모두 지우기



file.txt라는 이름의 파일을 남기고 나머지를 모두 지우기.


find . ! -name 'file.txt' -type f -exec rm -f {} +



디렉토리라면 type을 d로 바꾸고 rm 에 recursive 옵션을 추가하면 된다.


find . ! -name 'file.txt' -type d -exec rm -rf {} +

* 주의사항 

위의 디렉토리만 남기고 모두 지우기는 디렉토리 안의 파일도 모두 지움

디렉토리 수준에서만 적용하고 싶다면 아래처럼 사용한다


find . -maxdepth 1 ! -name 'file.txt' -type d -exec rm -rf {} +
depth제한을 둠으로써 하위 폴더로 내려가지 않고 현재 폴더에서만 적용.


여러 개의 파일을 남겨놓고 싶다면 아래처럼 추가 할 수 있다.


find . ! -name 'file.txt' ! -name 'file2.txt' -type f -exec rm -f {} +


Reference -

https://unix.stackexchange.com/questions/153862/remove-all-files-directories-except-for-one-file

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

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

linux 계정 정보 옮기기  (1) 2019.08.21
tar 디렉토리 지정해서 압축 풀기  (0) 2018.10.04
리눅스에서 프록시 설정하기  (0) 2018.07.26
neocomplcache vim plugin 설치하기  (0) 2018.07.12
samba 설정하기  (0) 2018.02.02
반응형

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

+ Recent posts

티스토리툴바