Fold change에 따른 DEG 분석(Python)

바이오마커

• 표현형(phenotype)의 표지/원인이 되는 유전자
• 질병의 진단, 치료 등에 활용

 

차별적 발현 유전자(DEG, Differentially Expressed Genes)

• RNA 데이터에서 차별적 발현 유전자(DEG) 검출 분석
  • 과발현 : 실험군(e.g. 폐암환자군)에서 이상적으로 양이 많이 발현 되어 있는 유전자
  • 저발현 : 실험군(e.g. 폐암환자군)에서 이상적으로 양이 적게 발현 되어 있는 유전자
• 차별적 발현 유전자는 꼭 바이오마커라고는 볼 수 없지만 바이오마커의 후보라고 볼 수 있음

 

DEG를 찾아내는 3가지 방법

• Fold Change(계산이 단순, 간단하게 분석하고 싶을 때, 옛날에 많이 썼음, 학술적이지 않다는 의견)
• P-value(Statistics, 통계학적 접근)
• Machine Learning(최근)

 

DEG by Fold Change

• 유전자 발현 데이터에서 실험/대조군의 평균 차이 비율(fold change)를 구함
• 보통 fold-change > 2, fold-change < 0.5 인 유전자를 DEG로 결정(상황에 따라 임의로 조정)

 

Fold Change

Fold change is a measure describing how much a quantity changes going from an initial to a final value.

• Normal Cell: 1, 2, 2, 2, 3 (평균: 2)
• Cancer Cell: 4, 4, 6, 6 (평균:5)

인 상황을 가정해보자. 

• Pseudo count 1을 도입하지 않으면 5/2 = 2.5
• Pseudo count 1을 도입하면 (5+1)/(2+1) = 2이다. (분모가 0이 되어 값이 무한대가 되는 것을 방지하기 위해 도입)

=> 2 이상이므로  DEG라고 판단

구해진 FC를 그냥 쓰기도 하지만 보통 log2를 씌워서 양수->증가, 음수->감소라고 쉽게 판단하기도 한다.

 

Lower-bound expression 제거

평균적으로 값이 너무 작을 때 FC가 아무리 크더라도 DEG로 뽑지 않는 경우가 있음.

Genes should be expressed > X at least Y% samples

관습적으로 1 이상이 발현된 유전자가 80% 이상이 되어야 함.

• Normal Cell: 0, 0, 0, 0, 2
• Cancer Cell: 0, 0, 2, 3, 0 

인 상황을 가정해보자. 전체 10개 중 30%의 sample에서만 발현이 되었다.

이는 lower bound(80%)를 넘지 못했기 때문에 FC의 값에 상관 없이 DEG로 뽑지 않음.

 

실습 파일 및 코드

http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/LUAD/20160128/gdac.broadinstitute.org_LUAD.Merge_rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.Level_3.2016012800.0.0.tar.gz

import pandas as pd
data=pd.read_csv("LUAD.txt",delimiter="\t",skiprows=lambda x: x == 1,index_col=0)

lst_labels = []
for col in data.columns:
    # if col = "TCGA-05-4244-01A-01R-1107-07", then sample_type = "01A"
    sample_type = col.split('-')[3] 
    if sample_type.startswith("01"): # primary tumor
        lst_labels.append('cancer')
    elif sample_type.startswith("02"): # reccurent solid tumor
        lst_labels.append('cancer')
    elif sample_type.startswith("10"): # blood derived normal
        lst_labels.append('normal')
    elif sample_type.startswith("11"): # solid tissue normal
        lst_labels.append('normal')
    else:
        lst_labels.append('normal')
        print("Warnning: sample type out of options")
        
        import numpy as np
lst_log2fold_change = []
lst_overone_ratio=[]

for i in range(data.shape[0]):
    # load i-th row (i.e, data of gene i)
    line= list(data.iloc[i,:])
    
    # For computing log2 fold change of gene i
    lst_cancer_vals=[]
    lst_normal_vals=[]
    for j in range(len(line)):
        if lst_labels[j] == 'cancer':
            lst_cancer_vals.append(line[j])
        elif lst_labels[j] == 'normal':
            lst_normal_vals.append(line[j])
    log2fold_change = np.log2((np.mean(lst_cancer_vals)+1) / (np.mean(lst_normal_vals)+1))
    lst_log2fold_change.append(log2fold_change)

    # For computing over one ration of gene i
    lst_overone=[]
    for j in range(len(line)):
        if line[j] > 1:
            lst_overone.append(1)
        else:
            lst_overone.append(0)
    overone_ratio = sum(lst_overone)/len(lst_overone)
    lst_overone_ratio.append(overone_ratio)
    lst_isDEG=[]
    
for i in range(data.shape[0]):
    log2fold_change = lst_log2fold_change[i]
    overone_ratio = lst_overone_ratio[i]
    if log2fold_change > np.log2(1.8) and overone_ratio > 0.8:
        isDEG = 1
    elif log2fold_change < -np.log2(1.8) and overone_ratio > 0.8:
        isDEG = -1
    else:
        isDEG = 0
    lst_isDEG.append(isDEG)
res_final = pd.DataFrame({'gene':data.index,'log2fold_change':lst_log2fold_change,'FCDEG':lst_isDEG})
res_final.to_csv("LUAD.foldchange", sep='\t', header=True, index=False)