import scanpy as sc
import harmonypy as hm
# 1. 데이터 로드
# 10개의 데이터 파일 경로 지정
file_paths = [
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606920_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606921_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606922_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606923_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606924_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606925_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606926_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606927_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606928_analysis/outs/filtered_feature_bc_matrix",
"/data/project/kim89/Cell_ranger/cellranger-8.0.1/SRR11606929_analysis/outs/filtered_feature_bc_matrix"
]
# 각 데이터를 Scanpy의 AnnData로 읽어들임
adatas = [sc.read_10x_mtx(path, var_names="gene_symbols", cache=True) for path in file_paths]
# donor 정보 추가
# donor ID 리스트
donor_ids = [
"donor_1", "donor_2", "donor_3",
"donor_4_1", "donor_5_1",
"donor_4_2", "donor_5_2",
"donor_6", "donor_7", "donor_8"
]
# 각 데이터에 donor ID를 추가
for adata, donor_id in zip(adatas, donor_ids):
adata.obs["batch"] = donor_id
# 데이터 확인
# for adata in adatas:
# print(adata.obs["batch"].unique())
# 2. 데이터 결합
# 공통 유전자 기준으로 데이터를 결합
# 모든 데이터를 결합
merged_adata = adatas[0].concatenate(
*adatas[1:], batch_key="batch", batch_categories=donor_ids
)
# 결합된 데이터의 batch 확인
# print(merged_adata.obs["batch"].unique())
print("Merged data:", merged_adata)
"""
Merged data: AnnData object with n_obs × n_vars = 37686 × 38606
obs: 'batch'
var: 'gene_ids', 'feature_types'
"""
# 3. QC (필터링 전 기본 QC 메트릭 계산)
# 미토콘드리아 유전자 필터링
merged_adata.var["mt"] = merged_adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(merged_adata, qc_vars=["mt"], inplace=True)
# before QC plot
sc.settings.figdir = "/data/project/kim89/SRR116069__plots"
sc.pl.violin(
merged_adata,
keys=["total_counts", "n_genes_by_counts", "pct_counts_mt"],
jitter=0.4,
multi_panel=True,
save="violin_plot_QC_before.png"
)
# QC 필터링 기준 설정 (Fan 기준 사용)
QC_filter = (
(merged_adata.obs["n_genes_by_counts"] > 200) & # 유전자 수 > 200
(merged_adata.obs["n_genes_by_counts"] < 2500) & # 유전자 수 < 2500
(merged_adata.obs["total_counts"] > 300) & # UMI > 300
(merged_adata.obs["total_counts"] < 15000) & # UMI < 15000
(merged_adata.obs["pct_counts_mt"] < 10) # 미토콘드리아 비율 < 10%
)
# 필터링 적용
filtered_adata = merged_adata[QC_filter, :]
print("Filtered data:", filtered_adata)
"""
Filtered data: View of AnnData object with n_obs × n_vars = 33062 × 38606
obs: 'batch', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'
"""
# after QC plot
# # Fan 기준
sc.pl.violin(
filtered_adata,
keys=["total_counts", "n_genes_by_counts", "pct_counts_mt"],
jitter=0.4,
multi_panel=True,
save="violin_plot_QC_after.png"
)
# 필터링된 데이터 저장 (.h5ad 형식으로 저장하여 이후 분석에 사용)
# filtered_adata.write_h5ad("filtered_adata.h5ad")
# 정규화 전 데이터의 총 카운트 분포 확인
# 33062개의 세포, 38606개의 유전자
# total_counts의 평균: 약 3045 # 분포가 세포마다 다르고, 최대값이 약 9037로 세포 간 차이가 큼.
print("정규화 전:")
print(filtered_adata.obs["total_counts"].describe())
"""
정규화 전:
count 33062.000000
mean 3047.830322
std 1660.008179
min 500.000000
25% 1942.000000
50% 2775.000000
75% 3787.000000
max 14487.000000
Name: total_counts, dtype: float64
"""
# 정규화 전 첫 번째 세포의 값
print("정규화 전 첫 번째 세포:")
print(merged_adata.X[0, :])
"""
정규화 전 첫 번째 세포:
<Compressed Sparse Row sparse matrix of dtype 'float32'
with 1290 stored elements and shape (1, 38606)>
Coords Values
(0, 58) 1.0
(0, 81) 1.0
...
"""
# 4. 정규화
# 총 발현값으로 정규화
sc.pp.normalize_total(filtered_adata, target_sum=1e4)
# 로그 변환
sc.pp.log1p(filtered_adata)
# 변수 유전자 선택
sc.pp.highly_variable_genes(filtered_adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
filtered_adata = filtered_adata[:, filtered_adata.var.highly_variable]
# 데이터 확인
# 정규화 후 총합 다시 계산
sc.pp.calculate_qc_metrics(filtered_adata, inplace=True)
# 정규화 후 데이터의 총 카운트 분포 확인
# 33062개의 세포, 3052개의 변수 유전자 (38606개의 유전자에서 줄어듦)
# total_counts의 평균: 약 272 # 모든 세포가 동일한 총합(target_sum)에 맞추어 정규화되었으며, 값의 범위가 좁아짐.
print("정규화 후:")
print(filtered_adata.obs["total_counts"].describe())
"""
정규화 후:
count 33062.000000
mean 272.048615
std 89.915779
min 40.955204
25% 214.527405
50% 265.780319
75% 319.822426
max 1039.758545
Name: total_counts, dtype: float64
"""
# 정규화 후 첫 번째 세포의 값
print("정규화 후 첫 번째 세포:")
print(filtered_adata.X[0, :])
"""
정규화 후 첫 번째 세포:
<Compressed Sparse Row sparse matrix of dtype 'float32'
with 196 stored elements and shape (1, 3052)>
Coords Values
(0, 8) 1.6668635606765747
(0, 26) 2.2608320713043213
"""
print(filtered_adata)
"""
AnnData object with n_obs × n_vars = 33062 × 3052
obs: 'batch', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'feature_types', 'mt', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
uns: 'log1p', 'hvg'
"""
# normalization plot
sc.pl.violin(
filtered_adata,
keys=["total_counts", "n_genes_by_counts", "pct_counts_mt"],
jitter=0.4,
multi_panel=True,
save="violin_plot_normalization.png"
)
# 5. 배치 효과 조정 (Harmony 사용), PCA, 클러스터링
import harmonypy as hm
# PCA 계산
sc.tl.pca(filtered_adata)
print(filtered_adata.obsm['X_pca'].shape) # (셀 수, PCA 차원 수) : (33062, 50)
# Harmony 실행
harmony_results = hm.run_harmony(
filtered_adata.obsm['X_pca'], # PCA 결과
filtered_adata.obs, # 관측값 데이터 (obs) (셀 수, 메타데이터 열 수)
'batch' # 배치 효과를 조정할 열 이름
)
# Harmony 결과 확인
print("Before transpose:", harmony_results.Z_corr.shape) # (50, 33062)
# Z_corr는 (PCA 차원, 셀 수) 형식으로 반환되므로 이를 Transpose(전치)하여 (셀 수, PCA 차원) 형식으로 맞춰야 합니다.
# Transpose하여 저장
filtered_adata.obsm['X_pca_harmony'] = harmony_results.Z_corr.T
print("After transpose, Harmony 결과 저장 완료:", filtered_adata.obsm['X_pca_harmony'].shape) # (33062, 50)
# Neighbors 계산
sc.pp.neighbors(filtered_adata, use_rep="X_pca_harmony")
print("Neighbors 계산 완료")
# UMAP 계산
sc.tl.umap(filtered_adata)
print("UMAP 계산 완료")
# Leiden 클러스터링 수행
sc.tl.leiden(filtered_adata, resolution=1.0) # resolution 값 조정 가능
print("Leiden 클러스터링 결과:", filtered_adata.obs['leiden'].value_counts())
"""
Leiden 클러스터링 결과: leiden
0 6462
1 4474
2 4293
3 3766
4 3062
5 2473
6 1813
7 1748
8 1532
9 997
10 974
11 886
12 174
13 165
14 133
15 66
16 44
Name: count, dtype: int64
"""
# UMAP 시각화
sc.pl.umap(filtered_adata, color=["batch", "leiden"], save="batch_effect_umap.png")
print("UMAP 시각화 완료")
# 데이터 저장
# filtered_adata.write_h5ad("merged_filtered_normalized_batch_corrected.h5ad")
print("데이터 저장 완료")
############# << 3. 탐색분석 >> ###########
# 0. DEGs
# 클러스터 간 차이를 확인하기 위해 주요 차별 발현 유전자(DEGs)를 식별
sc.tl.rank_genes_groups(filtered_adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(filtered_adata, n_genes=20, save="ranked_genes.png")
# 클러스터별 UMAP 시각화
sc.pl.umap(filtered_adata, color=["leiden", "PF4"], save="umap_gene_highlight_PF4.png") ##### 특정 유전자 NKG7의 cluster 존재여부
# 1. 기능 강화 분석 (Functional Enrichment Analysis)
# 생물학적 기능적 특성을 도출할 수 있습니다. 이를 위해 Gene Ontology (GO), KEGG Pathway 분석, 그리고 GSVA를 사용할 수 있습니다.
# (1) 주요 유전자 집합 준비
# 주요 유전자 추출
import pandas as pd
top_genes = pd.DataFrame(filtered_adata.uns['rank_genes_groups']['names']).head(50) # 상위 50개 유전자
# print(top_genes)
# 0 1 2 ... 14 15 16
# 0 LEF1 LYZ NKG7 ... GNLY CD74 HSP90AA1
# 1 CAMK4 CST3 GNLY ... NKG7 AFF3 NKG7
# [50 rows x 17 columns]
# 클러스터 0의 상위 유전자 추출
cluster_0_genes = filtered_adata.uns['rank_genes_groups']['names']['0'][:50] # 클러스터 0의 상위 50개 유전자
# print(cluster_0_genes)
# ['LEF1' 'CAMK4' 'TSHZ2' 'PRKCA' 'INPP4B' 'PDE3B' 'SERINC5' 'BCL11B' 'FHIT'
# 'NELL2' 'BACH2' 'TXK' 'BCL2' 'ANK3' 'ENSG00000271774' 'TC2N' 'TRABD2A'
# 'PLCL1' 'ENSG00000290067' 'MLLT3' 'ENSG00000249806' 'LNCRNA-IUR' 'SESN3'
# 'LINC01550' 'IGF1R' 'THEMIS' 'PATJ' 'PRKCA-AS1' 'ENSG00000253980' 'NR3C2'
# 'LINC02964' 'RETREG1' 'CSGALNACT1' 'APBA2' 'EDA' 'ABLIM1' 'ZC3H12B' 'AK5'
# 'ITGA6' 'RASGRF2' 'IMMP2L' 'ANKRD55' 'PRKN' 'LEF1-AS1' 'ICOS' 'KANK1'
# 'MDS2' 'ZEB1' 'ENSG00000286330' 'FAAH2']
# print(filtered_adata.uns['rank_genes_groups']['names']['0'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['1'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['2'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['3'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['4'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['5'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['6'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['7'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['8'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['9'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['10'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['11'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['12'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['13'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['14'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['15'][:10])
# print(filtered_adata.uns['rank_genes_groups']['names']['16'][:10])
# (2) GO 및 KEGG Pathway 분석
# gprofiler-official를 사용한 분석
from gprofiler import GProfiler
# GO 및 KEGG 분석
gp = GProfiler(return_dataframe=True)
# 클러스터 0의 상위 유전자 리스트로 변환
cluster_0_genes_list = list(set(cluster_0_genes.tolist())) # 중복 제거 포함
print(type(cluster_0_genes)) # <class 'numpy.ndarray'>
print(type(cluster_0_genes.tolist())) # <class 'list'>
# gProfiler 분석
go_results = gp.profile(organism='hsapiens', query=cluster_0_genes_list) ##### cluster 0에서의 GO
# print(go_results.head()) # 상위 몇 줄 출력
# source native name ... recall query parents
# 0 GO:BP GO:0048856 anatomical structure development ... 0.004238 query_1 [GO:0032502]
# 1 GO:BP GO:0032502 developmental process ... 0.004029 query_1 [GO:0008150]
# 2 GO:BP GO:0010646 regulation of cell communication ... 0.005518 query_1 [GO:0007154, GO:0050794]
# 3 GO:BP GO:0007275 multicellular organism development ... 0.004523 query_1 [GO:0032501, GO:0048856]
# 4 GO:BP GO:0009653 anatomical structure morphogenesis ... 0.005961 query_1 [GO:0032502, GO:0048856]
# [5 rows x 14 columns]
# print(go_results.columns) # 열 이름 확인
# Index(['source', 'native', 'name', 'p_value', 'significant', 'description',
# 'term_size', 'query_size', 'intersection_size', 'effective_domain_size',
# 'precision', 'recall', 'query', 'parents'],
# dtype='object')
print(go_results[['name', 'p_value', 'source']]) # 결과 확인
# 저장
go_results.to_csv("GO_KEGG_enrichment_results.csv", index=False)
import matplotlib.pyplot as plt
import seaborn as sns
# p-value를 기준으로 상위 10개 경로 시각화
top_results = go_results.nsmallest(10, 'p_value')
sns.barplot(y=top_results['native'], x=-np.log10(top_results['p_value']))
plt.xlabel('-log10(p-value)')
plt.ylabel('GO Terms / KEGG Pathways')
plt.show()
# from pyGSVA import gsva
# # Gene set 정의 (MSigDB 또는 KEGG)
# gene_sets = {
# "Apoptosis": ["TP53", "BCL2", "CASP3"], # 예시 유전자
# "Cell_Cycle": ["CDK1", "CDK2", "CCNB1"]
# }
# # GSVA 분석
# gsva_results = gsva(expression_data=filtered_adata.to_df(), gene_sets=gene_sets, method='gsva')
# print(gsva_results)
##########
# # 2. 전사인자 분석 (Transcription Factor Analysis)
# import scanpy as sc
# from pyscenic.pipeline import run_pyscenic_pipeline
# # 데이터 준비
# filtered_adata.raw = filtered_adata # raw data 저장 (pySCENIC은 raw counts 필요)
# filtered_adata.X = filtered_adata.raw.to_adata().X # raw counts 사용
# # pySCENIC 실행
# run_pyscenic_pipeline(
# adata=filtered_adata,
# output_prefix='pyscenic_results',
# organism='homo_sapiens',
# dbs=['cisTarget_human_v10_clust.genes_vs_motifs.rankings.db'],
# n_threads=4
# )
# # 결과 파일 읽기
# import pandas as pd
# regulons = pd.read_csv('pyscenic_results/regulons.csv')
# # 상위 전사인자 시각화
# sc.pl.dotplot(filtered_adata, var_names=regulons.head(20), groupby='leiden')
# # 3. 의사시간 분석 (Pseudotime Analysis)
# # 의사시간 계산
# sc.tl.dpt(filtered_adata, n_dcs=10) # DPT 계산
# sc.pl.umap(filtered_adata, color='dpt_pseudotime', save="pseudotime.png")
# # 4. 세포 간 상호작용 (Cell-Cell Interaction Analysis)
# import squidpy as sq
# # Spatial data 준비 (spatial data가 없을 경우 random assignment 가능)
# adata = sq.datasets.imc() # 예제 데이터
# # Ligand-Receptor 상호작용 계산
# sq.gr.ligrec(adata)
# # 시각화
# sq.pl.ligrec(adata, target_groups="leiden", source_groups="batch", dendrogram=True)
# # 5.세포 주기 분석 (Cell Cycle Analysis)
# # 세포 주기 점수 계산
# s_genes = ["MCM5", "PCNA", "TYMS", "FEN1", "MCM2", "MCM4", "RRM1", "UNG"]
# g2m_genes = ["CCNB1", "CCNB2", "CDC20", "CDK1", "PLK1", "TOP2A"]
# sc.tl.score_genes_cell_cycle(filtered_adata, s_genes=s_genes, g2m_genes=g2m_genes)
# # UMAP 시각화
# sc.pl.umap(filtered_adata, color=["phase", "S_score", "G2M_score"], save="cell_cycle.png")
