总结：
-
数据读取：
- 加载了多个数据文件，包括RNA-seq数据、等位基因表达数据、样本信息以及BAM文件的统计数据。
-
等位基因数据处理：
- 通过
separate()、spread()和rename()等函数对等位基因数据进行了长宽格式转换，确保每个样本的两个等位基因的计数值分别存储在pa1和pa2列中。 - 使用
replace_na()将缺失值替换为零。
- 通过
-
计算等位基因比例：
- 对每个样本计算了
pa1的比例，表示该等位基因在该样本中的相对丰度。
- 对每个样本计算了
-
合并与展示：
- 计算结果通过
inner_join()与样本元数据（如基因型、处理类型等）进行了合并。 - 最终打印出前120行的合并数据，帮助进一步分析样本的特征与基因表达情况。
- 计算结果通过
这个脚本的核心目的是对等位基因的表达比例进行计算并与样本元数据结合，最终得到每个样本的等位基因表达比例以及相关的元数据信息，方便后续的分析与研究。
第一部分：数据读取和初始化
source("functions.R")
t_cfg %>% print(n=50)yid = 'rn17c'
yid = 'rn18g'
diri = file.path(dird, 'raw', yid)
fi = file.path(diri, 'cpm.rds')
x = readRDS(fi)
fi = file.path(diri, 'ase.rds')
ase = readRDS(fi)
th = x$th %>% select(SampleID,Tissue,Genotype,Treatment)
fi = file.path(diri, 'ase0.tsv')
ase0 = read_tsv(fi)
fi = file.path(diri, 'bamstats.tsv')
bs = read_tsv(fi)
解析：
-
source("functions.R")
这一行加载了一个外部的R脚本文件functions.R，其中可能定义了一些用于数据处理的自定义函数。通过调用source()，你可以在当前脚本中使用该文件中的所有函数。 -
t_cfg %>% print(n=50)
t_cfg可能是一个配置对象或数据框，通过管道操作符%>%将其传递给print()函数，打印前50行数据。t_cfg可能包含了实验配置或元数据。 -
yid = 'rn17c'和yid = 'rn18g'
这两行代码似乎是设置了不同的实验或样本的ID。实际运行时，只有第二个yid = 'rn18g'会被有效地赋值（因为它覆盖了前一个）。 -
diri = file.path(dird, 'raw', yid)
创建了一个文件路径diri，指向存放原始数据的目录。dird可能是一个根目录变量，yid是前面设置的样本ID（rn18g），这条路径指向与样本ID相关的文件夹。 -
fi = file.path(diri, 'cpm.rds')
创建了指向cpm.rds文件的路径。这个文件存储了RNA-seq数据（可能是经过计算后的每百万计数值，CPM）。 -
x = readRDS(fi)
使用readRDS()读取cpm.rds文件，并将其存储在变量x中。x可能是一个包含RNA-seq数据的列表或数据框，其中包含了基因表达信息。 -
fi = file.path(diri, 'ase.rds')
创建了指向ase.rds文件的路径。ase可能是一个与RNA-seq数据相关的文件，包含了某些分析结果（例如，表达的等位基因信息）。 -
ase = readRDS(fi)
读取ase.rds文件并将其赋值给变量ase，它存储了基因组数据（可能与等位基因表达量相关）。 -
th = x$th %>% select(SampleID, Tissue, Genotype, Treatment)
从x中提取th数据框（通常包含样本元数据），并选择了四列（SampleID、Tissue、Genotype、Treatment）。这些列存储了与每个样本相关的元数据（样本ID、组织类型、基因型和处理类型）。 -
fi = file.path(diri, 'ase0.tsv')和ase0 = read_tsv(fi)
创建了指向ase0.tsv文件的路径并读取该文件，将其存储在ase0变量中。ase0可能包含了其他形式的基因组数据，格式为TSV。 -
fi = file.path(diri, 'bamstats.tsv')和bs = read_tsv(fi)
创建了指向bamstats.tsv文件的路径，并使用read_tsv()读取该文件，存储在bs变量中。bamstats.tsv可能包含与BAM文件相关的统计信息，如读取的数量、覆盖度等。
总结：
这一部分主要是加载了实验数据，包括RNA-seq的表达数据（cpm.rds），等位基因表达数据（ase.rds），样本元数据（x$th），以及额外的基因组数据（ase0.tsv、bamstats.tsv）。这些数据为后续的分析（如基因表达模式分析、样本信息处理等）提供了必要的输入。
很好！接下来是第二部分：对等位基因数据进行处理。
第二部分：对等位基因数据进行处理
ase1 = ase %>%separate(sid, c('sid', 'pa'), sep='[\\.]') %>%spread(pa, cnt) %>%rename(pa1 = `1`, pa2 = `2`) %>%replace_na(list(pa1=0, pa2=0))
解析：
-
ase %>% separate(sid, c('sid', 'pa'), sep='[\\.]')- 这里，
ase是之前从ase.rds读取的等位基因数据。 separate()函数用于将sid列按分隔符（[\\.]，即“.”）拆分成两列：sid和pa。sid可能是样本的ID，而pa代表的是等位基因的不同部分（例如，基因的两个等位基因可能会通过这个列标记）。
- 这里，
-
spread(pa, cnt)spread()函数用于将数据从“长格式”转换为“宽格式”。- 这里，将
pa列（表示等位基因的部分）转换为列名，并将cnt（计数值）分配到相应的列。结果就是，等位基因的计数（cnt）将被分配到pa1和pa2两列，分别表示两个等位基因的计数。
-
rename(pa1 =1, pa2 =2)rename()函数将列名1和2分别重命名为pa1和pa2，这使得数据更具可读性和直观性。这里的pa1和pa2代表两个等位基因的计数值。
-
replace_na(list(pa1=0, pa2=0))replace_na()函数用于替换缺失值（NA）。- 在这里，如果
pa1或pa2存在缺失值，则将其替换为0，表示这些等位基因没有被检测到或计数为0。
总结：
这一部分主要是对 ase 数据进行处理，将原始的等位基因数据从长格式转换为宽格式，并处理缺失值。经过处理后，ase1 数据框将包含每个样本的两个等位基因（pa1 和 pa2）的计数。
好的，接下来是第三部分：计算等位基因比例并与样本元数据合并。
第三部分：计算等位基因比例并与样本元数据合并
ase1 %>% group_by(sid) %>%summarise(prop.B = sum(pa1) / sum(pa1 + pa2)) %>%ungroup() %>%inner_join(th, by = c('sid' = 'SampleID')) %>%print(n=120)
解析：
-
ase1 %>% group_by(sid)- 这行代码将
ase1数据框按sid（样本ID）进行分组。每个组代表一个样本，之后会对每个样本的等位基因计数进行聚合计算。
- 这行代码将
-
summarise(prop.B = sum(pa1) / sum(pa1 + pa2))summarise()函数用于生成每个组（每个样本）的汇总统计信息。- 在这里，通过计算
pa1和pa2（代表两个等位基因的计数）的和，求出pa1占总计数（pa1 + pa2）的比例，结果存储在新列prop.B中。prop.B表示每个样本中，等位基因pa1的比例。
-
ungroup()ungroup()解除之前的分组操作，这样后续的操作（如inner_join()）不再基于分组进行。这样可以防止在后续操作中意外地保持分组状态。
-
inner_join(th, by = c('sid' = 'SampleID'))- 使用
inner_join()将计算得到的每个样本的等位基因比例（prop.B）与之前加载的样本元数据（th）合并。 - 合并的条件是
sid（来自ase1）与SampleID（来自th）相匹配。这样，每个样本的等位基因比例会和样本的其他信息（如组织类型、基因型、处理类型等）一起出现在合并后的数据框中。
- 使用
-
print(n=120)- 打印合并后的数据框的前120行，以便查看计算结果和合并后的样本信息。通过这种方式，你可以检查每个样本的等位基因比例以及与其相关的元数据。
总结：
这一部分的代码计算了每个样本中 pa1 的比例（即 prop.B），然后将这些信息与样本的元数据合并，最终打印出前120行合并后的数据。合并后，每个样本的等位基因比例与其他信息（如基因型、处理类型等）一起展示，便于后续分析和理解样本的特征。
你清楚这一部分了吗？如果理解了，我们可以继续。