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

Python

Gene Expression Omnibus

scRNA-seq

Bioinformatics

Author

Taeyoon Kim

Published

November 9, 2024

Modified

June 30, 2025

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 anndata as ad
import numpy as np
import pandas as pd
from scipy.io import mmread
from scipy.sparse import csr_matrix


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


def load_tsv_file(
    file_path: str, header: int | None = None, names: list[str] | None = 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) -> ad.AnnData | None:
    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: ad.AnnData | None = process_sample(data_dir, sample_name)
        if adata is not None:
            adatas.append(adata)

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

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

adatas[0].obs.head()

adatas[0].var.head()

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

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()

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

adata_merged.obs.shape

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