我有一个表,我试图根据一列中基于另一列中的变量的特定均值子集对其进行规范化。理想情况下,我的代码应该将特定应变变量(如2987)的coverage_depth列中的所有数据除以同一列的一个子集的平均值(仅针对chr列中的SAG1的覆盖深度,仅针对应变列中的2987 )
我已经找到了一种很长的方法来做这件事,但我真的希望有人能把它变成一个循环,这样我就不需要在计算完均值后手动输入它们。
我的表看起来像这样:
B1 1073 320 2987
B1 1074 324 2987
B1 1075 330 2987
SAG1 955 31 2987
SAG1 956 30 2987
SAG1 957 29 2987
SAG1 958 29 2987
BTub 446 57 2987
BTub 452 59 2987
B1 1707 53 GRE_MIG
B1 1708 56 GRE_MIG
18S 1099 242 GRE_MIG
18S 1100 242 GRE_MIG
SAG1 888 7 GRE_MIG
SAG1 889 7 GRE_MIG
SAG1 890 7 GRE_MIG
首先,我在我的表中加载:
reads<-read.table("3133_all.CNV.txt", sep = "\t", header = F)
colnames(reads)<-c("chr", "position", "coverage_depth", "strains"
然后,我调用plyr来计算所有chr和coverage_depth列的组合的平均值。
library(plyr)
coverage_summary<-ddply(reads, c("chr", "strains"), summarise, mean = mean(coverage_depth))
write.csv(format(coverage_summary, scientific=FALSE), file = "CNV_mean_07.27.16.csv", row.names = F)
这给了我一个更长的版本:
chr strains mean
1 18S 2987 2.052802e+03
20 18S GRE_MIG 2.674536e+01
126 B1 GRE_MIG 6.503342e+01
213 SAG1 2987 3.422057e+01
232 SAG1 GRE_MIG 5.863501e+00
我想出了如何对一个菌株的所有coverage_depth进行归一化,这是我在chr SAG1中从该菌株获得的平均值,我手动输入了如下内容:
NormalizeSAG1<-function(coverage_depth, strains){
if (strains %in% c("2987")) {
coverage_depth<-coverage_depth/3.42
} else if (strains %in% c("GRE_MIG")) {
coverage_depth<-coverage_depth/5.86
} else { coverage.norm<-coverage_depth
}}
reads$SAG1_normalized<-mapply(NormalizeSAG1, reads$coverage_depth, reads$strains)
问题是我有53个不同的菌株,我想根据它们在chr列中单个SAG1的平均值来归一化它们。似乎一个for循环可以做到这一点,但是我不知道如何在没有大量ifelse语句的情况下正确地将我的数据子集以进行规范化。
发布于 2016-08-09 23:46:47
尝试以下操作:
reads <- merge(reads, coverage_summary)
reads <- mutate(reads, normalized = coverage_depth / mean)
基本上,这应该将您的摘要列连接回原始数据,之后,创建一个规范化列应该是微不足道的。这也避免了必须创建一个自定义函数来处理53个不同的可能值。
https://stackoverflow.com/questions/38862035
复制