LDSC不仅仅可以用于计算
以及检验GWAS结果的
inflation是否是因为混杂因素（confounding）的影响，还可以有S-LDSC（partition haritability， 分区遗传力），LDSC-SEG（遗传力组织细胞富集）
在本节，让我们来讲讲LDSC-SEG这个方法，实质上也是利用了partition haritability的思想。
一、S-LDSC的理论及方法
S-LDSC的思想是通过不同的基因组注释（不同功能区域）来估计在基因组上的分布。其拟合回归模型为
其中为SNP
的卡方统计量期望值，
为样本量，
为混杂因素，
是SNP
中的LD评分，
为待估计的回归系数。
-
- 若
- 若
- 若
- 若
通过回归，我们可以得到每个注释类别中的SNP的富集，及哪些注释类别对于表型变异更加重要。
二、S-LDSC的脚本实现
1、转换格式
python munge_sumstats.py \
--sumstats body_BMIz.sumstats.gz \
--merge-alleles w_hm3.snplist \
--out BMI
2、运行S-LDSC
python ldsc/ldsc.py \
--h2 BMI.sumstats.gz \
--ref-ld-chr baselineLD. \
--frqfile-chr 1000G.EUR.QC. \
--w-ld-chr weights.hm3_noMHC. \
--overlap-annot \
--print-coefficients \
--print-delete-vals \
--out BMI.baselineLD
上面脚本结束后，我们就可以得到S-LDSC的结果，查看在哪些注释类别中更富集。
三、LDSC-SEG组织细胞富集
继S-LDSC的思想，研究团队又将partition haritability拓展到了转录组的组织细胞层面。
1、首先，计算某个数据来源（比如这里的GTEx）表达矩阵 per gene 的
t-statistic。
2、对某个组织或细胞的t-statistic进行排序，比如这里的cortex组织。
3、根据t-statistic排序，取top 10%的gene作为该组织活细胞specific expression的gene。
4、对top 10%的gene上下100-KB窗口，利用LDSC baseline model计算其LD score，并利用S-LDSC进行组织细胞富集。
LDSC-SEG.png
！！！注意：如果得到的基因表达矩阵为标准化后的相对表达矩阵， 则不用计算t-statistic，就用该值进行排序选基因即可。
The available gene expression values already quantify relative expression for a tissue or cell type rather than absolute expression for a single sampleand so we used these values in place of our t-statistics. (Hilary K. Finucane et al. Nat Genet)
四、LDSC-SEG代码实现
1、转换格式
python munge_sumstats.py \
--sumstats body_BMIz.sumstats.gz \
--merge-alleles w_hm3.snplist \
--out BMI
2、运行LDSC-SEG
ldsc.py \
--h2-cts BMI.sumstats.gz \
--ref-ld-chr 1000G_EUR_Phase3_baseline/baseline. \
--out BMI_${cts_name} \
--ref-ld-chr-cts $cts_name.ldcts \
--w-ld-chr weights_hm3_no_hla/weights.
五、如何生成自己的LDSC-SEG组织细胞注释文件
因为LDSC只提供了一些常见组织的注释文件，那么如果我们研究的是一些如耳朵、鼻子等组织，就得自己找数据生成这些注释文件。有时，有些在人类中难以得到，则需要用小鼠这些的代替，不过需要先找同源基因。\
我这里就用小鼠的scRNA数据进行实例,前提是你已经有了每个细胞类型的基因表达矩阵。
1、biomaRt找同源基因
# if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# BiocManager::install("biomaRt")
library(dplyr)
library(data.table)
library(biomaRt)
listEnsembl()
# connect Ensembl
mouse <- useEnsembl(biomart = "genes", dataset = "mmusculus_gene_ensembl")
human <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
## 构建转化函数
homologs <- getBM(attributes = c("chromosome_name","ensembl_gene_id", "external_gene_name",
"mmusculus_homolog_ensembl_gene", "mmusculus_homolog_orthology_type"),
filters = "with_mmusculus_homolog",
values = TRUE,
mart = human)
# head(homologs)
homologs_1to1 <- subset(homologs, mmusculus_homolog_orthology_type == "ortholog_one2one")
# rm MT X Y chromosome genes
homGene_1to1_chr1_22 <- homologs_1to1 %>%
filter(!chromosome_name %in% c("MT", "X", "Y")) %>%
rename(human_ensembl_gene = ensembl_gene_id)
fwrite(homGene_1to1_chr1_22, file="HomGene1to1_mouse2human.txt", sep="\t", row.names=FALSE)
2、得到per cell type top 10%的基因list
# 同时生成control.GeneSet文件，也就是所有的基因
library(dplyr)
library(data.table)
setwd("/public/share/wchirdzhq2022/Wulab_share/LDSC/Mouse_cochleae/")
scGene = fread("Sc_P8_12_20.txt")
homologs_1to1 = fread("HomGene1to1_mouse2human.txt")[,2:3,with=FALSE]
colnames(homologs_1to1) = c("human_gene", "gene")
merge_df = merge(scGene, homologs_1to1, by="gene")
merge_df = subset(merge_df, select = -c(gene))
colnames(merge_df) = gsub("[[:space:]/-]", "_", colnames(merge_df))
colnames(merge_df) = gsub("'", "_", colnames(merge_df))
merge_df = merge_df %>% select(human_gene, everything())
control = merge_df[, 1]
fwrite(control, file="Sc_P8_12_20_control.GeneSet", sep="\t", col.names=FALSE)
lapply(2:ncol(merge_df), function(i){
temp_dt = data.frame(merge_df[[1]], merge_df[[i]])
# colnames(df1) = names(merge_df)[c(37, i)]
sorted_dt = temp_dt[order(temp_dt[[2]], decreasing = TRUE), ]
top_10 = sorted_dt[1:floor(0.1*nrow(sorted_dt)), ]
top_10_gene = top_10[, 1, drop=FALSE]
cell_type = colnames(merge_df)[i]
fwrite(top_10_gene, file=paste0("Sc_P8_12_20_", cell_type, ".GeneSet"), sep="\t", col.names=FALSE)
})
3、make_annot.py处理得到的gene list
虽然LDSC说make_annot.py默认提供了ENSG.coord.txt，但是其实这个文件找不到，得自己生成。
由于这个文件需要的是GENE CHR START END四列信息，而SMR中提供的glist-hg19文件正好满足，所以用这个文件稍微转化一下
# generate ENSG.coord.txt
library(biomaRt)
library(data.table)
human <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
human_gene <- getBM(attributes = c("chromosome_name","ensembl_gene_id", "external_gene_name"),
filters = "", values = TRUE, mart = human)
colnames(human_gene) = c("CHR", "ENSG_ID", "GENE")
glist_hg19 <- fread("/public/home/shilulu/script/plot_smr/glist-hg19")
colnames(glist_hg19) = c("CHR", "START", "END", "GENE")
merge_df <- merge(glist_hg19, human_gene, by="GENE")
glist_hg19_ENSG <- merge_df[, c("ENSG_ID", "CHR.y", "START", "END")]
colnames(glist_hg19_ENSG) = c("GENE", "CHR", "START", "END")
fwrite(glist_hg19_ENSG, file="ENSG_coord.txt", sep="\t")
然后，运行make_annot.py，其中的control为你的整个表达矩阵中所有的基因哦，control在运行LDSC-SEG的时候会用到
# using 100K widow as Hilary K. Finucane et al. Nat Genet
# run make_annot.py 为control 和 每个cell type生成
conda activate ldsc
bfile="/public/share/wchirdzhq2022/Wulab_share/LDSC/1000G_EUR_Phase3_plink/1000G.EUR.QC"
ldsc="/public/home/shilulu/software/ldsc"
ls *.GeneSet | while read id; do
annot=$(basename -- "${id}" ".GeneSet")
for chr in {1..22}; do
cmd="python ${ldsc}/make_annot.py \
--gene-set-file ${annot}.GeneSet \
--gene-coord-file ${ldsc}/ENSG_coord.txt \
--bimfile ${bfile}.${chr}.bim \
--windowsize 100000 \
--annot-file ${annot}.${chr}.annot.gz"
done
done
4、计算per cell type的LD score
awk '{if ($1!="SNP") {print $1} }' ${ldsc_dir}/w_hm3.snplist > ${ldsc_dir}/listHM3.txt
ls *.GeneSet | while read id; do
annot=$(basename -- "${id}" ".GeneSet")
for chr in {1..22}; do
python ${ldsc}/ldsc.py --l2 \
--bfile ${bfile}.${chr} \
--print-snps ${ldsc_dir}/listHM3.txt \
--ld-wind-cm 1 \
--annot ${annot}.{TASK_ID}.annot.gz \
--thin-annot \
--out ${annot}.${chr}
done
done
5、生成一个list，第一列为cell type name，第二列为每个cell type文件名前缀加“,”加上面说的cotrol的文件名前缀
# generate .ldcts file for my sc data
ls ${dir}/*GeneSet | while read id; do
tissue=$(basename -- ${id} | sed 's/^Sc_P8_12_20_//;s/.GeneSet$//')
cell=${dir}/$(basename -- ${id} "GeneSet")
control=${dir}/"Sc_P8_12_20_control."
echo "${tissue} ${cell},${control}" >> Mouse_cochleae_gene_expr.ldcts
done
# 注意最后把cell type为cotrol的删除，这个脚本有待完善
至此，我们就生成了自己的组织或细胞类型的注释文件，就可以进行LDSC-SEG分析了
ldsc.py --h2-cts BMI.sumstats.gz \
--ref-ld-chr ${REF_LD_CHR} \
--out BMI_cts \
--ref-ld-chr-cts Mouse_cochleae_gene_expr.ldcts \
--w-ld-chr ${W_LD_CHR}
