如何用wcAggregateCases解决多通道序列分析的序列数量超限问题？

问题背景 #

我有54399个样本，包含HOM和HOS两个通道，计划开展多通道序列分析，数据示例如下：

我使用了以下代码：

HOM.seq<-seqdef(df[, 2:5])
HOS.seq<-seqdef(df[, 6:9])
channels<-list(HOM.seq, HOS.seq)
MDdist<-seqMD(channels, method="OM", sm=list("TRATE", "TRATE"), what="diss")

但运行后收到警告：52322 unique sequences exceeds max allowed of 46340

我的问题：

• 如何使用wcAggregateCases减少唯一序列的数量？尽管这52322个序列似乎已经是从54399个样本聚合而来。
• 是否可以在将HOM和HOS放入通道列表前，分别使用wcAggregateCases进行聚合？我曾单独对HOM和HOS使用wcAggregateCases，聚合后的序列数量分别约为10000和7000。

解决方案 #

1. 针对多通道序列整体聚合 #

你看到的52322个唯一序列是两个通道组合后的唯一序列数，不是单个通道的聚合结果。要对多通道序列整体用wcAggregateCases，需要先合并两个通道的序列对象，再执行聚合：

# 创建多通道序列对象
multi_seq <- seqbind(HOM.seq, HOS.seq, varname=c("HOM", "HOS"))
# 聚合多通道序列，得到聚合后的序列及对应权重
agg_multi <- wcAggregateCases(multi_seq)
agg_seq <- agg_multi$aggseq
weights <- agg_multi$weights

# 基于聚合后的序列计算距离
MDdist <- seqMD(list(agg_seq[,1:4], agg_seq[,5:8]), method="OM", sm=list("TRATE", "TRATE"), what="diss", weights=weights)

2. 先单独聚合单通道再组合的可行性 #

完全可以先对单个通道分别聚合，再组合计算距离，这种方式能灵活控制每个通道的聚合粒度：

# 聚合HOM通道
agg_HOM <- wcAggregateCases(HOM.seq)
HOM_agg <- agg_HOM$aggseq
HOM_w <- agg_HOM$weights

# 聚合HOS通道
agg_HOS <- wcAggregateCases(HOS.seq)
HOS_agg <- agg_HOS$aggseq
HOS_w <- agg_HOS$weights

# 组合成多通道列表，计算组合权重（单个通道权重相乘）
channels_agg <- list(HOM_agg, HOS_agg)
combined_weights <- HOM_w * HOS_w

# 传入权重计算距离
MDdist_agg <- seqMD(channels_agg, method="OM", sm=list("TRATE", "TRATE"), what="diss", weights=combined_weights)

注意：单独聚合后组合时，权重需要用两个通道的权重相乘——因为每个聚合后的HOM序列对应HOM_w个原样本，每个聚合后的HOS序列对应HOS_w个原样本，组合后的多通道序列对应的样本数是两者的乘积。

额外优化技巧 #

如果想进一步压缩序列数量，可以通过minfreq参数过滤低频序列，比如只保留出现次数≥5的序列：

agg_HOM <- wcAggregateCases(HOM.seq, minfreq=5)

另外，聚合后计算距离时必须传入权重，否则会丢失原样本的数量信息，导致分析结果偏差。

内容的提问来源于stack exchange，提问作者Fanny0000