多物种比较转录组¶

# === OrthoFinder直系同源鉴定 ===

# 准备各物种的蛋白序列文件
mkdir proteomes
cp human.pep.fa proteomes/Human.fa
cp mouse.pep.fa proteomes/Mouse.fa
cp zebrafish.pep.fa proteomes/Zebrafish.fa
cp chicken.pep.fa proteomes/Chicken.fa

# 运行OrthoFinder
orthofinder -f proteomes/ -t 16 -a 8

# 关键输出文件：
# Results/Orthogroups/Orthogroups.tsv          - 所有orthogroup
# Results/Orthogroups/Orthogroups_SingleCopyOrthologues.txt - 1:1直系同源
# Results/Orthologues/Human__v__Mouse.tsv      - 成对同源关系
# Results/Gene_Trees/                          - 每个OG的基因树

# 提取1:1同源对（最适合跨物种比较）
# 1:1同源：每个物种只有一个基因的orthogroup
cat Results/Orthogroups/Orthogroups_SingleCopyOrthologues.txt | wc -l

# 解析OrthoFinder结果，构建跨物种基因映射
import pandas as pd

# 读取orthogroups
og = pd.read_csv("Orthogroups.tsv", sep="\t", index_col=0)

# 提取1:1:1:1同源基因（所有物种各一个）
one_to_one = og.dropna()  # 没有缺失物种
one_to_one = one_to_one[one_to_one.map(lambda x: ',' not in str(x)).all(axis=1)]
# 排除含多个paralogs的行
# 注意：pandas >= 2.1 中 applymap 已弃用，改用 map

print(f"1:1:1:1 orthologs: {len(one_to_one)}")
# 通常在脊椎动物间有~10000-13000个1:1同源

# 创建基因ID映射表
gene_mapping = pd.DataFrame({
    'Human': one_to_one['Human'].str.strip(),
    'Mouse': one_to_one['Mouse'].str.strip(),
    'Zebrafish': one_to_one['Zebrafish'].str.strip(),
    'Chicken': one_to_one['Chicken'].str.strip()
})
gene_mapping.to_csv("ortho_gene_mapping.csv", index=True)

# === 跨物种表达量准备与归一化 ===
library(dplyr)

# 读取各物种的TPM矩阵
human_tpm <- read.csv("human_tpm.csv", row.names = 1)
mouse_tpm <- read.csv("mouse_tpm.csv", row.names = 1)
zebrafish_tpm <- read.csv("zebrafish_tpm.csv", row.names = 1)

# 读取基因映射
mapping <- read.csv("ortho_gene_mapping.csv", row.names = 1)

# 提取1:1同源基因的表达量
human_ortho <- human_tpm[mapping$Human, , drop = FALSE]
mouse_ortho <- mouse_tpm[mapping$Mouse, , drop = FALSE]
zebrafish_ortho <- zebrafish_tpm[mapping$Zebrafish, , drop = FALSE]

# 统一行名为orthogroup ID
rownames(human_ortho) <- rownames(mapping)
rownames(mouse_ortho) <- rownames(mapping)
rownames(zebrafish_ortho) <- rownames(mapping)

# log2转换
human_log <- log2(human_ortho + 1)
mouse_log <- log2(mouse_ortho + 1)
zebrafish_log <- log2(zebrafish_ortho + 1)

# 方法1：Quantile归一化跨物种
library(preprocessCore)
# 合并后统一quantile normalize
all_expr <- cbind(human_log, mouse_log, zebrafish_log)
all_qn <- normalize.quantiles(as.matrix(all_expr))
colnames(all_qn) <- colnames(all_expr)
rownames(all_qn) <- rownames(all_expr)

# 方法2：中位数对齐
# 每个物种的每个样本减去中位数，加上全局中位数
global_median <- median(c(as.matrix(human_log), as.matrix(mouse_log)))
human_aligned <- sweep(human_log, 2, apply(human_log, 2, median) - global_median)
mouse_aligned <- sweep(mouse_log, 2, apply(mouse_log, 2, median) - global_median)

# === 表达保守性分析 ===

# 假设有多组织数据：brain, liver, heart, kidney
# 计算同源组织间的表达相关性

# 人与鼠brain表达相关
cor_brain <- cor(human_log[, "brain"], mouse_log[, "brain"],
                  method = "spearman")
cat(sprintf("Human-Mouse brain Spearman r = %.3f\n", cor_brain))

# 全组织表达谱相关矩阵
tissues <- c("brain", "liver", "heart", "kidney")
cor_matrix <- matrix(NA, length(tissues), length(tissues))
rownames(cor_matrix) <- paste0("Human_", tissues)
colnames(cor_matrix) <- paste0("Mouse_", tissues)

for (i in seq_along(tissues)) {
  for (j in seq_along(tissues)) {
    cor_matrix[i, j] <- cor(human_log[, tissues[i]], mouse_log[, tissues[j]],
                             method = "spearman", use = "complete.obs")
  }
}

# 热图可视化
library(pheatmap)
pheatmap(cor_matrix, display_numbers = TRUE,
         main = "Cross-species tissue expression correlation")
# 期望：对角线最高（同源组织最相似）

# 每个基因的跨物种保守性评分
gene_conservation <- sapply(rownames(human_log), function(gene) {
  h_expr <- as.numeric(human_log[gene, tissues])
  m_expr <- as.numeric(mouse_log[gene, tissues])
  cor(h_expr, m_expr, method = "spearman")
})

# 保守基因 vs 快进化基因
highly_conserved <- names(sort(gene_conservation, decreasing = TRUE)[1:500])
rapidly_evolving <- names(sort(gene_conservation)[1:500])

# 功能富集
library(clusterProfiler)
ego_conserved <- enrichGO(gene = mapping$Human[rownames(mapping) %in% highly_conserved],
                           OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP")
ego_divergent <- enrichGO(gene = mapping$Human[rownames(mapping) %in% rapidly_evolving],
                           OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP")

# === 表达进化速率 ===

# 方法1：表达散度（Expression Divergence）
# 使用1-Spearman r作为表达距离
expr_divergence <- 1 - gene_conservation

# 方法2：EVE模型（Expression Variance and Evolution）
# EVE区分种内变异（noise）和种间进化（真实变化）
# 类似于McDonald-Kreitman检验的表达版本

# 方法3：Ornstein-Uhlenbeck模型
# 使用phylogenetic方法建模表达进化
# 检验是否有表达最优值的约束（稳定化选择）

# 表达散度 vs 序列散度的关系
# 读取Ka/Ks值
kaks <- read.csv("pairwise_kaks.csv", row.names = 1)
merged_evol <- data.frame(
  gene = rownames(kaks),
  ka_ks = kaks$Ka_Ks,
  expr_div = expr_divergence[rownames(kaks)]
)

# 散点图
library(ggplot2)
ggplot(merged_evol, aes(x = ka_ks, y = expr_div)) +
  geom_point(alpha = 0.3, size = 0.5) +
  geom_smooth(method = "lm") +
  labs(x = "Ka/Ks (sequence divergence)",
       y = "Expression divergence (1 - Spearman r)") +
  theme_minimal()

# 通常正相关但弱——表达和序列进化部分解耦
cor.test(merged_evol$ka_ks, merged_evol$expr_div, method = "spearman")

# === PAML/codeml 正选择检测 ===

# 前提：需要多物种CDS比对和物种树

# 1. 准备CDS比对（codon-aware alignment）
# 使用PRANK或MACSE做保留密码子结构的比对
prank -d=orthogroup.cds.fa -t=species_tree.nwk -o=aligned -codon

# 2. 转换为PAML格式
# codeml需要phylip格式的比对文件

# 3. 编写codeml控制文件（位点模型: M7 vs M8）
cat > codeml_site.ctl << 'EOF'
      seqfile = aligned.phy       * 序列比对文件
     treefile = species_tree.nwk  * 物种树
      outfile = site_results.txt  * 输出文件
        noisy = 3
      verbose = 1
      runmode = 0                 * 0=用户提供树
      seqtype = 1                 * 1=密码子
    CodonFreq = 2                 * 2=F3x4模型
        model = 0                 * 0=所有分支同一omega
      NSsites = 7 8              * 7=M7(beta), 8=M8(beta+正选择)
        icode = 0                 * 0=标准遗传密码
    fix_omega = 0                 * 0=估计omega
        omega = 0.5               * 初始omega值
EOF

# 运行codeml
codeml codeml_site.ctl

# 4. 分支-位点模型（检测特定谱系的正选择）
cat > codeml_branchsite.ctl << 'EOF'
      seqfile = aligned.phy
     treefile = species_tree_marked.nwk  * 前景分支用#1标记
      outfile = branchsite_results.txt
      runmode = 0
      seqtype = 1
    CodonFreq = 2
        model = 2                 * 2=分支-位点模型
      NSsites = 2                 * 2=正选择
    fix_omega = 0                 * 备选假说：omega自由估计
        omega = 1.5
EOF
# 零假说：fix_omega=1, omega=1（不允许正选择）
# 用似然比检验(LRT)比较两个模型：2*(lnL1-lnL0)服从卡方分布(df=1)
# BEB(Bayes Empirical Bayes)方法识别具体受正选择的位点

codeml codeml_branchsite.ctl

# === 解析codeml结果 ===
# M7 vs M8比较
# 从输出文件提取lnL值
# LRT statistic = 2 * (lnL_M8 - lnL_M7)
# 自由度 = 2（M8比M7多2个参数）
# p = pchisq(LRT, df=2, lower.tail=FALSE)

# BEB结果中 posterior probability > 0.95 的位点为正选择位点

# === 物种特异性表达变化检测 ===

# 方法：在系统发育背景下识别某物种偏离祖先状态的基因

# 简单方法：Z-score
# 对每个基因，计算人类表达在所有物种中的Z-score
# Z > 2: 人类特异性高表达
# Z < -2: 人类特异性低表达

species_mean_expr <- data.frame(
  Human = rowMeans(human_log),
  Mouse = rowMeans(mouse_log),
  Zebrafish = rowMeans(zebrafish_log)
)

# 对每个基因计算人类的Z-score
species_mean_expr$human_zscore <- scale(species_mean_expr$Human)

# 人类特异性高表达基因（相对于其他物种）
for (i in 1:nrow(species_mean_expr)) {
  vals <- as.numeric(species_mean_expr[i, 1:3])
  species_mean_expr$human_zscore[i] <- (vals[1] - mean(vals)) / sd(vals)
}

human_specific_up <- rownames(species_mean_expr)[species_mean_expr$human_zscore > 2]
human_specific_down <- rownames(species_mean_expr)[species_mean_expr$human_zscore < -2]

cat(sprintf("Human-specific upregulated: %d\n", length(human_specific_up)))
cat(sprintf("Human-specific downregulated: %d\n", length(human_specific_down)))

# === 表达距离树 ===

# 合并所有物种-组织的平均表达谱
all_profiles <- cbind(
  human_log[, tissues],
  mouse_log[, tissues],
  zebrafish_log[, tissues]
)
colnames(all_profiles) <- c(
  paste0("Human_", tissues),
  paste0("Mouse_", tissues),
  paste0("Zebrafish_", tissues)
)

# 计算距离矩阵
dist_matrix <- as.dist(1 - cor(all_profiles, method = "spearman"))

# 层次聚类
hc <- hclust(dist_matrix, method = "average")
plot(hc, main = "Expression-based clustering of species-tissues",
     xlab = "", sub = "")

# 观察聚类模式：
# 如果Brain_Human先与Brain_Mouse聚，说明脑的表达谱跨物种保守
# 如果Human_Brain先与Human_Liver聚，说明物种效应>组织效应
# 通常结果：高度保守组织（脑、睾丸）跨物种聚类，其他按物种聚类

# === Paralog表达分化 ===

# 提取含paralogs的orthogroups
# 从OrthoFinder结果中找人类有两个基因的orthogroup
og_all <- read.table("Orthogroups.tsv", sep = "\t", header = TRUE, row.names = 1)
human_paralogs <- og_all[grepl(",", og_all$Human), ]

# 对每对paralog计算表达相关性
paralog_cors <- numeric()
for (i in 1:nrow(human_paralogs)) {
  genes <- strsplit(human_paralogs$Human[i], ", ")[[1]]
  if (length(genes) == 2 && all(genes %in% rownames(human_tpm))) {
    r <- cor(as.numeric(human_tpm[genes[1], ]),
             as.numeric(human_tpm[genes[2], ]),
             method = "spearman")
    paralog_cors <- c(paralog_cors, r)
  }
}

# 与1:1同源基因的保守性比较
hist(paralog_cors, breaks = 50, main = "Paralog Expression Correlation",
     xlab = "Spearman r", col = "lightblue")
abline(v = median(gene_conservation), col = "red", lty = 2,
       lwd = 2)  # 同源基因保守性作为参考

# 表达高度分化的paralog（neo-functionalization的候选）

# 同源基因鉴定
orthofinder -f proteomes/ -t 16 -a 8

# 各物种RNA-seq标准分析
hisat2 -x species_index -1 R1.fq.gz -2 R2.fq.gz | samtools sort -o aligned.bam
stringtie aligned.bam -G species.gtf -e -B -o quant.gtf
# 或salmon定量
salmon quant -i species_index -l A -1 R1.fq.gz -2 R2.fq.gz -o salmon_out