GSE에서 sc-RNAseq 데이터 가져오기

GSE 데이터베이스는 고처리량 유전자 발현 데이터, 특히 유전자 발현 옴니버스(Gene Expression Omnibus) 데이터의 저장소입니다. 이는 국립생물공학정보센터(NCBI)에서 관리하는 무료 온라인 데이터베이스로, 연구자와 과학자들이 유전자 발현 데이터를 공유하고 접근하는 데 널리 사용됩니다. 이 데이터베이스는 RNA 시퀀싱 실험을 포함한 다양한 출처의 방대한 데이터셋 컬렉션을 포함하고 있고 연구자들이 다운로드하여 추가 분석과 연구 목적으로 활용할 수 있습니다. 이는 생물정보학과 유전체학 연구에 매우 귀중한 자원입니다.

1 예시 데이터

예시 데이터로 GSE231559를 불러오는 방법을 다뤄 보겠습니다. GSE231559는 2023년 9월 25일에 공개된 데이터셋으로, 대장암(CRC)의 단일 세포 수준에서의 유전자 발현 프로파일을 다루고 있습니다. 이 연구는 MD Anderson 암센터의 CRC Moon Shot 프로젝트의 일환으로 수행되었으며, 단일 세포 RNA 시퀀싱(scRNA-seq) 기술을 사용하여 대장암의 다양성을 분석했습니다. 이 데이터셋은 Oncogenic KRAS drives lipo–fibrogenesis to promote angiogenesis and colon cancer progression 논문에 사용되었습니다.

2 원시 데이터 살펴보기

GSE 데이터베이스에서 파일을 다운로드해서 압축을 풀어보면 GSE231559가 아래와 같이 여러개의 파일로 되어 있다는 것을 확인 할 수 있습니다.

GSM7290760_SC10_21N_barcodes.tsv.gz  GSM7290773_SC10_5_barcodes.tsv.gz
GSM7290760_SC10_21N_features.tsv.gz  GSM7290773_SC10_5_features.tsv.gz
GSM7290760_SC10_21N_matrix.mtx.gz    GSM7290773_SC10_5_matrix.mtx.gz
GSM7290761_SC10_22T_barcodes.tsv.gz  GSM7290774_SC10_7_barcodes.tsv.gz
GSM7290761_SC10_22T_features.tsv.gz  GSM7290774_SC10_7_features.tsv.gz
GSM7290761_SC10_22T_matrix.mtx.gz    GSM7290774_SC10_7_matrix.mtx.gz
GSM7290762_SC10_23N_barcodes.tsv.gz  GSM7290775_SC10_8_barcodes.tsv.gz
GSM7290762_SC10_23N_features.tsv.gz  GSM7290775_SC10_8_features.tsv.gz
GSM7290762_SC10_23N_matrix.mtx.gz    GSM7290775_SC10_8_matrix.mtx.gz
GSM7290763_SC10_24T_barcodes.tsv.gz  GSM7290776_SC143_17_barcodes.tsv.gz
GSM7290763_SC10_24T_features.tsv.gz  GSM7290776_SC143_17_features.tsv.gz
GSM7290763_SC10_24T_matrix.mtx.gz    GSM7290776_SC143_17_matrix.mtx.gz
GSM7290764_SC10_25N_barcodes.tsv.gz  GSM7290777_SC143_7_barcodes.tsv.gz
GSM7290764_SC10_25N_features.tsv.gz  GSM7290777_SC143_7_features.tsv.gz
GSM7290764_SC10_25N_matrix.mtx.gz    GSM7290777_SC143_7_matrix.mtx.gz
GSM7290765_SC10_26N_barcodes.tsv.gz  GSM7290778_SC173_1_barcodes.tsv.gz
GSM7290765_SC10_26N_features.tsv.gz  GSM7290778_SC173_1_features.tsv.gz
GSM7290765_SC10_26N_matrix.mtx.gz    GSM7290778_SC173_1_matrix.mtx.gz
GSM7290766_SC10_27N_barcodes.tsv.gz  GSM7290779_SC173_2_barcodes.tsv.gz
GSM7290766_SC10_27N_features.tsv.gz  GSM7290779_SC173_2_features.tsv.gz
GSM7290766_SC10_27N_matrix.mtx.gz    GSM7290779_SC173_2_matrix.mtx.gz
GSM7290767_SC10_28T_barcodes.tsv.gz  GSM7290780_SC216_1_barcodes.tsv.gz
GSM7290767_SC10_28T_features.tsv.gz  GSM7290780_SC216_1_features.tsv.gz
GSM7290767_SC10_28T_matrix.mtx.gz    GSM7290780_SC216_1_matrix.mtx.gz
GSM7290768_SC10_29N_barcodes.tsv.gz  GSM7290781_SC216_2_barcodes.tsv.gz
GSM7290768_SC10_29N_features.tsv.gz  GSM7290781_SC216_2_features.tsv.gz
GSM7290768_SC10_29N_matrix.mtx.gz    GSM7290781_SC216_2_matrix.mtx.gz
GSM7290769_SC10_30T_barcodes.tsv.gz  GSM7290782_SC216_3_barcodes.tsv.gz
GSM7290769_SC10_30T_features.tsv.gz  GSM7290782_SC216_3_features.tsv.gz
GSM7290769_SC10_30T_matrix.mtx.gz    GSM7290782_SC216_3_matrix.mtx.gz
GSM7290770_SC10_35N_barcodes.tsv.gz  GSM7290783_SC216_5_barcodes.tsv.gz
GSM7290770_SC10_35N_features.tsv.gz  GSM7290783_SC216_5_features.tsv.gz
GSM7290770_SC10_35N_matrix.mtx.gz    GSM7290783_SC216_5_matrix.mtx.gz
GSM7290771_SC10_37N_barcodes.tsv.gz  GSM7290784_SC216_6_barcodes.tsv.gz
GSM7290771_SC10_37N_features.tsv.gz  GSM7290784_SC216_6_features.tsv.gz
GSM7290771_SC10_37N_matrix.mtx.gz    GSM7290784_SC216_6_matrix.mtx.gz
GSM7290772_SC10_38T_barcodes.tsv.gz  GSM7290785_SC216_7_barcodes.tsv.gz
GSM7290772_SC10_38T_features.tsv.gz  GSM7290785_SC216_7_features.tsv.gz
GSM7290772_SC10_38T_matrix.mtx.gz    GSM7290785_SC216_7_matrix.mtx.gz

골치가 벌써 아파지네요. 그래도 파일명에 일정한 패턴이 있기에 아래 파이썬 코드를 사용해 갹갹의 adata 객체를 만들어서 합치면 될 것 같습니다.

import os
import scanpy as sc
import pandas as pd
import numpy as np
import anndata as ad
from scipy.io import mmread
from scipy.sparse import csr_matrix
from typing import List, Dict, Any, Optional
import matplotlib.pyplot as plt

def load_mtx_file(file_path: str) -> csr_matrix:
    return mmread(file_path).transpose().tocsr()

def load_tsv_file(file_path: str, header: Optional[int] = None, names: Optional[List[str]] = None) -> pd.DataFrame:
    return pd.read_csv(file_path, sep="\t", header=header, names=names)

def create_anndata(X: csr_matrix, obs_names: List[str], var_names: np.ndarray, var_gene_ids: np.ndarray, code_name: str) -> ad.AnnData:
    adata = ad.AnnData(X=X)
    adata.obs_names = obs_names
    adata.var_names = var_names
    adata.var["gene_ids"] = var_gene_ids
    adata.obs['Code_name'] = code_name
    return adata

def process_sample(data_dir: str, sample_name: str) -> Optional[ad.AnnData]:
    matrix_file = os.path.join(data_dir, f"{sample_name}_matrix.mtx.gz")
    features_file = os.path.join(data_dir, f"{sample_name}_features.tsv.gz")
    barcodes_file = os.path.join(data_dir, f"{sample_name}_barcodes.tsv.gz")
    code_name = sample_name.split("_")[0]

    if not all(os.path.exists(f) for f in [matrix_file, features_file, barcodes_file]):
        print(f"Missing files for {sample_name}")
        return None

    try:
        X = load_mtx_file(matrix_file)
        gene_names = load_tsv_file(features_file, names=["gene_ids", "gene_symbols", "Gene_expression"])
        barcodes = load_tsv_file(barcodes_file, names=["barcode"])

        obs_names = [f"{code_name}_{b}" for b in barcodes["barcode"]]
        var_names = gene_names["gene_symbols"].values
        var_gene_ids = gene_names["gene_ids"].values

        adata = create_anndata(X, obs_names, var_names, var_gene_ids, code_name)
        print(f"Successfully read {sample_name}")
        return adata
    except Exception as e:
        print(f"Error reading {sample_name}: {str(e)}")
        return None

data_dir: str = "../input/GSE231559"
adatas: List[ad.AnnData] = []

for file in os.listdir(data_dir):
    if file.endswith("_matrix.mtx.gz"):
        sample_name: str = file.replace("_matrix.mtx.gz", "")
        adata: Optional[ad.AnnData] = process_sample(data_dir, sample_name)
        if adata is not None:
            adatas.append(adata)

print(f"Total number of processed samples: {len(adatas)}")

Successfully read GSM7290779_SC173_2
Successfully read GSM7290775_SC10_8
Successfully read GSM7290761_SC10_22T
Successfully read GSM7290770_SC10_35N
Successfully read GSM7290768_SC10_29N
Successfully read GSM7290777_SC143_7
Successfully read GSM7290776_SC143_17
Successfully read GSM7290766_SC10_27N
Successfully read GSM7290774_SC10_7
Successfully read GSM7290763_SC10_24T
Successfully read GSM7290784_SC216_6
Successfully read GSM7290767_SC10_28T
Successfully read GSM7290764_SC10_25N
Successfully read GSM7290782_SC216_3
Successfully read GSM7290780_SC216_1
Successfully read GSM7290783_SC216_5
Successfully read GSM7290769_SC10_30T
Successfully read GSM7290773_SC10_5
Successfully read GSM7290771_SC10_37N
Successfully read GSM7290772_SC10_38T
Successfully read GSM7290765_SC10_26N
Successfully read GSM7290760_SC10_21N
Successfully read GSM7290762_SC10_23N
Successfully read GSM7290781_SC216_2
Successfully read GSM7290785_SC216_7
Successfully read GSM7290778_SC173_1
Total number of processed samples: 26

성공적으로 데이터를 불어왔는지 확인하기 첫번째 객체의 데이터를 확인해봅니다.

adatas[0].obs.head()

	Code_name
GSM7290779_AAACCCAAGCTAATCC-1	GSM7290779
GSM7290779_AAACGAAAGTCTAGCT-1	GSM7290779
GSM7290779_AAAGAACAGCTCAGAG-1	GSM7290779
GSM7290779_AAAGAACCATGAATAG-1	GSM7290779
GSM7290779_AAAGGATAGAGGCTGT-1	GSM7290779

adatas[0].var.head()

	gene_ids
MIR1302-2HG	ENSG00000243485
FAM138A	ENSG00000237613
OR4F5	ENSG00000186092
AL627309.1	ENSG00000238009
AL627309.3	ENSG00000239945

문제없이 불러온 것 같습니다.

3 여러 데이터를 하나로 합치기

이제 여러 데이터를 하나의 객체로 합쳐보겠습니다. 추가로 메타 데이터를 찾아서 맵핑을 통해 환자 정보와 질병 정보를 추가합니다.

def make_var_names_unique(adatas: List[ad.AnnData]) -> List[ad.AnnData]:
    for adata in adatas:
        adata.var_names_make_unique()
    return adatas

def merge_anndata(adatas: List[ad.AnnData]) -> ad.AnnData:
    adata_merged = ad.concat(adatas, join='outer', axis=0, index_unique="-")
    all_var = pd.concat([adata.var for adata in adatas], axis=0, join="outer")
    all_var = all_var.loc[~all_var.index.duplicated(keep='first')]
    adata_merged.var = all_var.loc[adata_merged.var_names]
    adata_merged.obs['Code_name'] = np.concatenate([np.repeat(adata.obs['Code_name'].iloc[0], adata.n_obs) for adata in adatas])
    return adata_merged

def create_sample_dict() -> Dict[str, str]:
    return {
        "GSM7290760": "L1_N", "GSM7290761": "L1_T", "GSM7290762": "C1_N",
        "GSM7290763": "C1_T", "GSM7290764": "L2_N", "GSM7290765": "L3_N",
        "GSM7290766": "L4_N", "GSM7290767": "L4_T", "GSM7290768": "C2_N",
        "GSM7290769": "C2_T", "GSM7290770": "L5_N", "GSM7290771": "C3_N",
        "GSM7290772": "C3_T", "GSM7290773": "C4_T", "GSM7290774": "C5_T",
        "GSM7290775": "L6_T", "GSM7290776": "L7_N", "GSM7290777": "C6_T",
        "GSM7290778": "L8_T", "GSM7290779": "L8_T", "GSM7290780": "L9_N",
        "GSM7290781": "L9_T", "GSM7290782": "L10_T", "GSM7290783": "L11_T",
        "GSM7290784": "L11_N", "GSM7290785": "L12_T"
    }

def map_code_to_sample(code: str, sample_dict: Dict[str, str]) -> str:
    return sample_dict.get(code, None)

def create_metadata() -> pd.DataFrame:
    metadata = {
        'Patient': ['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'L1', 'L4', 'L6', 'L8', 'L9', 'L11'],
        'Location': ['Colon', 'Colon', 'Colon', 'Colon', 'Colon', 'Colon', 'Liver', 'Liver', 'Liver', 'Liver', 'Liver', 'Liver'],
        'Gender': ['M', 'M', 'F', 'M', 'M', 'F', 'F', 'M', 'F', 'M', 'M', 'M'],
        'Current_age': [47, 39, 65, 85, 72, 62, 64, 70, 68, 54, 56, 60],
        'Stage_at_collection': ['IV', 'IV', 'IV', 'II', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV'],
        'MS_status': ['MSS', 'MSS', 'MSS', 'MSI', 'MSS', 'MSS', 'MSS', 'MSS', 'MSS', 'MSS', 'MSS', 'MSS']
    }
    meta_df = pd.DataFrame(metadata)
    meta_df.set_index('Patient', inplace=True)
    return meta_df

def extract_patient_id(sample: str) -> str:
    return sample.split('_')[0]

def determine_tissue(sample: str) -> str:
    parts = sample.split('_')
    if parts[-1].endswith('T'):
        return 'tumor'
    elif parts[-1].endswith('N'):
        return 'normal'
    else:
        return 'unknown'

def add_metadata_to_anndata(adata: ad.AnnData, meta_df: pd.DataFrame) -> ad.AnnData:
    adata.obs['Patient'] = adata.obs['Sample'].apply(extract_patient_id)
    for column in meta_df.columns:
        adata.obs[column] = adata.obs['Patient'].map(meta_df[column])
    adata.obs['Tissue_type'] = adata.obs['Sample'].apply(determine_tissue)
    return adata

# Main execution
adatas = make_var_names_unique(adatas)
adata_merged = merge_anndata(adatas)

sample_dict = create_sample_dict()
adata_merged.obs['Sample'] = adata_merged.obs['Code_name'].apply(lambda x: map_code_to_sample(x, sample_dict))

meta_df = create_metadata()
adata_merged = add_metadata_to_anndata(adata_merged, meta_df)

# Check the result
adata_merged.obs.head()

	Code_name	Sample	Patient	Location	Gender	Current_age	Stage_at_collection	MS_status	Tissue_type
GSM7290779_AAACCCAAGCTAATCC-1-0	GSM7290779	L8_T	L8	Liver	M	54.0	IV	MSS	tumor
GSM7290779_AAACGAAAGTCTAGCT-1-0	GSM7290779	L8_T	L8	Liver	M	54.0	IV	MSS	tumor
GSM7290779_AAAGAACAGCTCAGAG-1-0	GSM7290779	L8_T	L8	Liver	M	54.0	IV	MSS	tumor
GSM7290779_AAAGAACCATGAATAG-1-0	GSM7290779	L8_T	L8	Liver	M	54.0	IV	MSS	tumor
GSM7290779_AAAGGATAGAGGCTGT-1-0	GSM7290779	L8_T	L8	Liver	M	54.0	IV	MSS	tumor

위의 출력을 보니 데이터가 정상적으로 합쳐진 것 같습니다. 데이터의 크기도 확인해보죠.

adata_merged.obs.shape

(39992, 15)

4 데이터 저장하기

이제 데이터 분석에 사용하기 위해 결과를 h5ad 파일로 저장합니다.

# 결과 저장
adata_merged.write('../output/241015_GSE231559_merged.h5ad', compression="gzip")

5 마치며

이번 글에서 우리는 GSE231559 데이터셋을 불러오고 처리하는 과정을 살펴보았습니다. 이러한 과정은 생물정보학 연구에서 매우 중요한 단계이지만, 몇 가지 도전적인 측면이 있습니다. 특히 수동으로 데이터를 합치는 과정과 메타데이터를 추가하는 부분은 모든 데이터셋마다 통일되지 않아 작업자가 직접 작성해야 한다는 점이 번거로울 수 있습니다. 각 데이터셋의 고유한 특성과 구조로 인해 이 과정을 완전히 자동화하기는 어렵기 때문입니다.

그럼에도 불구하고, 이러한 과정을 직접 수행해보는 것은 데이터의 구조와 특성을 깊이 이해하는 데 큰 도움이 됩니다. 또한, 이를 통해 데이터 처리 능력을 향상시킬 수 있으며, 향후 유사한 작업을 더 효율적으로 수행할 수 있게 될 것입니다.

이 글을 통해 독자 여러분이 GSE231559와 같은 데이터셋을 직접 다루어보고, 생물정보학 데이터 처리의 실제적인 측면을 경험해보시기를 바랍니다. 이러한 경험은 여러분의 연구나 프로젝트에 큰 도움이 될 것입니다.

6 참고

• 이 데이터를 포함하고 있는 다른 CRC scRNA-seq 데이터 모음: https://onlinelibrary.wiley.com/doi/10.1002/tox.24157