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"로 채우고 비어있지 않다면 값을 유지하고 지나간다.

Argument받아서 처리하기

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

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

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

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

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

argument가 0개 일때의 조건문.

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 {} +

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

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

16진수 RGB코드 알아내는법

화면 내에서 특정 색이 가진 RGB코드를 쉽게 알아내려면 EST soft에서 제공하는 알씨 프로그램이 가장 편하다. (https://www.altools.co.kr/download/alsee.aspx)

알씨 프로그램에서 [보기]->[픽셀정보] 를 클릭하면 마우스의 위치를 따라가면서 RGB값을 표시해준다. 0x는 16진법이라는 표시니 제외하고 뒤에 숫자+알파벳 조합이 색의 정보이다.

DESeq2에서 heatmap, PCA, MA, volcano plot 그리기

분석을 위한 DataSet이 필요하다. 예시에서는 여러 cell line에 몇 가지 drug를 각기 다른 dose로 처리했을 때의 결과를 분석하고자 한다.

DEG분석은 dose에만 초점이 맞춰지도록 디자인 하였다.

res변수는 결과가 adjust p-value 낮은순서대로 sorting하여 저장하였다.

Cancer cell line 정보 받기

https://cancer.sanger.ac.uk/cell_lines

Catalogue Of Somatic Mutations In Cancer (COSMIC) 홈페이지에서 cancer cell line에 대한 정보를 받는 법을 알아보려고 한다. 특히나 특정 cancer cell line으로 실험을 하였을 때 결과가 제대로 나왔는지 확인하기 위해 해당 cell line에 존재하는 variant를 검사해 볼 수 있다.

가장 먼저 홈페이지에서 genome version을 설정해주어야 한다. 이 후의 자료들은 다 위의 genome version을 기반으로 제공된다.

옆에 검색창에 "HeLa"를 검색해보자.

비슷한 이름의 다른 정보가 하나도 없기 때문에 Samples에 하나 나온다. 나머지 정보는 뻔한 것들이니 넘어가고 뒷부분에 존재하는 variants를 보면 HeLa가 가지고 있는 variants가 보인다.

Filters를 사용하면 더 specific한 변이들만 따로 고를 수 있다.

Reference - 

https://cancer.sanger.ac.uk/cell_lines

Optical duplicate와 Library duplicates

NGS 데이터 생산시 reads의 숫자를 증폭시키기 위해 PCR cycle을 돌리게 되는데 분석 과정에서 과도하게 증폭된 reads를 제거하기위해 remove duplicate과정을 거치게 된다. 

picard의 markduplicate 모듈을 사용했을 때 나온 결과의 예시이다. Optical Duplicates와 Not Optical Duplicates가 존재한다.

데이터에서 관측되는 duplicate는 두 종류이다.

optical duplicates

- read를 읽을 때 single cluster를 인접한 두 개의 cluster라고 인식했을 때 생긴다. 두 sequence는 굉장히 유사한 sequence를 가지게 된다. 일반적으로 N개의 염기가 같은 경우 optical duplicates로 분류하는데 N은 read 전체 일 수도 있고 약 50정도를 사용할 수도 있다. alignment를 하지 않아도 찾아낼 수 있다.

library duplicates

- PCR duplicates라고도 말하며 라이브러리 제작 준비 과정에서 생겨날 수 있다. 독립적인 스팟에서 일어나기 때문에 optical duplicates와는 달리 NGS장비 내부에서 인접한 cluster에서 발생하지 않지만 서열이 굉장히 유사하다는 점이 optical duplicates와의 차이이다. genome에 맵핑 이후에 찾아낼 수 있다.

FastQC에서 보이는 duplicates를 어떻게 해석해야 할까?

- FastQC에서 보이는 duplicate level이 sequencing이 잘 되었다 안되었다를 의미하는 것은 아니다. 단지 데이터가 일반적인 특성에 따라가는지 여부만을 확인시켜 줄 뿐이다. 

(선 하나 하나가 하나의 sample을 보여주고 있는 것이다.)

- FastQC에서의 duplicate level은 perfect match만을 보여주는 것이기 때문에 sequencing error나 전체적인 quality score가 낮은 경우, 실제로는 duplication이 많지만 그렇지 않게 보일 수 있음을 주의해야 한다.

Reference -

https://wiki.bits.vib.be/index.php/Q%26A_added_during_the_intro_to_NGS_data_analysis

http://proteo.me.uk/2011/05/interpreting-the-duplicate-sequence-plot-in-fastqc/