liftover하기

liftover란?

- 다른 genome에 맞게 결과 파일의 버전을 바꾸는 것을 말한다. 

- 사람의 genome은 hg19, GRCh37, GRCh38 등 여러 개가 존재하며 이에 따라 각 유전자의 위치도 조금씩 차이가 난다. 새롭게 genome을 만들 때마다 모든 정보를 새로 작성하는 것 보다 기존의 정보에서 달라진 위치만 수정하는 것이 비용이나 시간상으로 효율적일 것이다.

- 권장하지는 않지만 다른 종의 genome 간의 비교도 가능하다. 단 여기서는 서열 간의 차이가 크게 나기 때문에 정확도가 낮아 손실되는 정보가 있을 수 있음에 유의해야 한다.

CrossMap

liftover를 지원하는 프로그램은 여러 개 있지만, 여기에선 CrossMap을 소개하고자 한다.

CrossMap은 SAM/BAM, Wiggle/BigWig, BED, GFF/GTF, VCF 등 다양한 포맷의 파일을 지원하며 특히나 python module이기 때문에 설치가 매우 간단하다.

아래와 같이 입력하면 설치가 완료된다.

Python 2.7 이하에서만 작동하는 모듈이다.

dependency가 있지만, 특별히 어려운 모듈은 없어서 pip에서 알아서 설치해주니 크게 신경 쓸 필요는 없다.

설치가 제대로 되었으면 chain 파일이 필요하다.

chain 파일이란 두 genome 간의 변화된 부분이 작성된 파일이다. CrossMap 홈페이지에서 사람과 쥐의 chain 파일을 제공하고 있으며 USCS genome browser에 가면 종간의 chain 파일도 내려받을 수 있다.

CrossMap 홈페이지 : http://crossmap.sourceforge.net/

프로그램 사용법은 홈페이지에 자세하게 나와 있으니 요약하기만 하겠다.

input file이 hg18이고 output file이 hg19라면 chain file은 hg18Tohg19를 넣으면 된다.

Reference -

http://crossmap.sourceforge.net/

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 이하이다.

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

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

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

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

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

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

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

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

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

pheatmap으로 heatmap그리기

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

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

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

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

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

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

Reference -

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

R에서 Dataframe 합치기

main이 같은 dataframe끼리 합친다.

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

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

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

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