玉米中的元基因调控网络突出了功能上相关的调控相互作用。\functions.R

1. 加载必要的R包

require(devtools)
load_all("~/git/rmaize")
require(ape)
require(ggtree)
require(ggforce)
require(Rtsne)

load_all("~/git/rmaize")：加载本地路径中的rmaize包，通常是你自己或团队开发的一个包，这里假设是一个与玉米相关的分析包。

2. 定义工作路径

dirp = '~/projects/rnaseq'
dird = glue('{dirp}/data')
#dirc = '/scratch.global/zhoux379/rnaseq'
#t_cfg = read_gspread_master(lib='rnaseq')

说明：

dirp = '~/projects/rnaseq'：定义RNA-Seq项目的根目录路径，指向本地的rnaseq项目文件夹。
dird = glue('{dirp}/data')：使用glue包来动态拼接路径，dird代表数据存储目录。
#dirc = '/scratch.global/zhoux379/rnaseq'：注释掉的行，可能是备用路径。
#t_cfg = read_gspread_master(lib='rnaseq')：注释掉的行，可能是从Google Sheets读取配置信息，配置RNA-Seq分析所需的文件和参数。

这一部分的目的是设定一些路径，以便后续操作中能轻松访问数据和相关文件。

3. 定义`read_rnaseq`函数

read_rnaseq <- function(yid) {res = rnaseq_cpm(yid)th = res$th; tm = res$tmth = th %>% replace_na(list(Tissue='', Genotype='B73', Treatment='')) %>%mutate(Tissue=as.character(Tissue)) %>%mutate(Genotype=as.character(Genotype)) %>%mutate(Treatment=as.character(Treatment))...
}

对 th 数据进行缺失值替换和数据转换

replace_na：这是 tidyr 包中的函数，用于替换缺失值。这里对 th 数据框的三个列进行了缺失值替换：
- Tissue 列的缺失值被替换为 ''（空字符串）。
- Genotype 列的缺失值被替换为 'B73'（这是一个常见的玉米品种名，可能是默认的基因型）。
- Treatment 列的缺失值被替换为 ''（空字符串）。
mutate：这是 dplyr 包中的函数，用于改变数据框中的列或创建新的列。这里使用 mutate 函数将 Tissue、Genotype 和 Treatment 列的类型转换为字符型（as.character）。这些列可能原本是因子型或其他类型，转换为字符型后，方便后续的处理和分析。

总结：

功能：该函数的核心目的是处理RNA-Seq数据中的 th 数据框，将缺失值进行合理替换，并确保某些列的格式一致（字符型）。这通常是数据预处理的一部分，为后续分析做准备。
输入：函数的输入 yid 可能是一个标识符，用于指定特定的RNA-Seq数据或实验条件。
输出：修改后的 th 数据框（包含了样本的组织、基因型和处理方式等信息）。

4. 定义`read_multiqc_trimmomatic`函数

read_multiqc_trimmomatic <- function(fi, paired = T) {ti = read_tsv(fi)types = c("surviving", "forward_only_surviving", "reverse_only_surviving", "dropped")if (paired == F) {nd = ti %>% mutate(nd = input_reads - surviving - dropped) %>%group_by(1) %>% summarise(nd = sum(nd)) %>% pull(nd)stopifnot(nd == 0)to = ti %>% mutate(SampleID = Sample, total = input_reads, surviving_f=0, surviving_r=0)} else if(paired == T | paired == 'both') {ti2 = ti %>% separate(Sample, c("SampleID",'suf'), sep="_", fill='right', extra='merge') %>% select(-suf) %>%mutate(surviving_f = forward_only_surviving,surviving_r = reverse_only_surviving)if(paired == 'both') ti2 = ti2 %>% replace_na(list('input_reads'=0, 'input_read_pairs'=0,'surviving_f'=0, 'surviving_r'=0)) %>%mutate(input_read_pairs = ifelse(input_read_pairs == 0, input_reads, input_read_pairs))nd = ti2 %>% mutate(nd = input_read_pairs - surviving - surviving_f - surviving_r - dropped) %>%group_by(1) %>% summarise(nd = sum(nd)) %>% pull(nd)stopifnot(nd == 0)to = ti2 %>% mutate(total = input_read_pairs)} else {stop(sprintf("unsupported option: %s", paired))}to %>%select(SampleID, total, surviving, surviving_f, surviving_r, dropped)
}

说明：

这个函数用于读取和处理MultiQC报告中的Trimmomatic数据。Trimmomatic是一个常用的高通量测序数据的质量控制工具。该函数的输入是fi（MultiQC报告文件路径）和paired（一个布尔值，指示是否是成对数据）。

读取数据：
- ti = read_tsv(fi)：使用read_tsv函数读取以制表符分隔的文件，返回一个数据框ti。
根据是否成对读取数据：
- 如果paired == F（即单端数据），会计算丢失的读取数（nd），并将surviving_f和surviving_r置为0。
- 如果paired == T（即成对数据），则使用separate(Sample, ...)分割样本列，分别提取前向和反向的存活读取数。
处理丢失读取数（nd）：
- 根据读取的数据计算nd，即总读取数减去存活读取数和丢弃读取数。如果nd不为0，则会停止并抛出错误。
- 如果paired == 'both'，则还会处理input_read_pairs和surviving_f、surviving_r。
输出：
- 最终函数返回一个数据框to，包含样本ID、总读取数、存活读取数、前向存活读取数、反向存活读取数和丢弃读取数。

该函数主要用于从质量控制文件中读取、处理和整理RNA-Seq数据的质量控制信息，特别是针对配对（或单端）数据的读数情况进行处理。它提供了对不同类型数据的灵活处理，并输出了包含样本、读数总数、存活读数等信息的简明数据框。

5. 定义`read_multiqc_star`函数

read_multiqc_star <- function(fi, paired = T) {ti = read_tsv(fi)if(paired == F) {ti2 = ti %>% mutate(SampleID = Sample) %>% select(-Sample)} else {ti2 = ti %>% separate(Sample, c("SampleID",'suf'), sep="_", fill='right', extra='merge') %>% select(-suf) }ti2 = ti2 %>%transmute(SampleID = SampleID, total = total_reads,uniquely_mapped = uniquely_mapped,multimapped = multimapped + multimapped_toomany,unmapped = unmapped_mismatches + unmapped_tooshort + unmapped_other,nd = total - uniquely_mapped - multimapped - unmapped) stopifnot(sum(ti2$nd) < 1000)ti2 = ti2 %>% group_by(SampleID) %>%summarise(uniquely_mapped = sum(uniquely_mapped),multimapped = sum(multimapped),unmapped = sum(unmapped))types = c("uniquely_mapped", "multimapped", "unmapped")to = ti2 %>% select(SampleID, uniquely_mapped, multimapped, unmapped)to
}

说明：

该函数处理MultiQC报告中的STAR比对数据。STAR是一个用于RNA-Seq数据的快速比对工具。函数处理STAR比对的输出数据并返回相关的统计信息。

读取数据：
- ti = read_tsv(fi)：读取MultiQC报告中的STAR比对统计数据。
处理成对和单端数据：
- 如果paired == F，则处理单端数据。
- 如果paired == T，则分离Sample列（假设是成对的样本名）。
计算不同的比对情况：
- transmute(SampleID = SampleID, total = total_reads, ...)：计算总读取数、唯一比对数、多重比对数、未比对数等。
- nd = total - uniquely_mapped - multimapped - unmapped：计算丢失的读取数（nd）。
数据汇总：
- 使用group_by(SampleID)按样本ID分组，汇总每个样本的比对统计信息。
- 最后返回一个包含每个样本唯一比对数、多重比对数和未比对数的数据框。

这个函数的主要目的是处理STAR比对工具的输出，并汇总每个样本的比对统计信息。

6. 定义`read_multiqc_bwa`函数

read_multiqc_bwa <- function(fi, paired = T) {ti = read_tsv(fi)if(paired == F) {ti2 = ti %>% mutate(SampleID = Sample) %>% select(-Sample)} else {ti2 = ti %>% separate(Sample, c("SampleID",'suf'), sep="_", fill='right', extra='merge') %>% select(-suf) }ti2 = ti2 %>%transmute(SampleID = SampleID, total = total_reads,uniquely_mapped = uniquely_mapped,multimapped = multimapped + multimapped_toomany,unmapped = unmapped_other,nd = total - uniquely_mapped - multimapped - unmapped) stopifnot(sum(ti2$nd) < 1000)ti2 = ti2 %>% group_by(SampleID) %>%summarise(uniquely_mapped = sum(uniquely_mapped),multimapped = sum(multimapped),unmapped = sum(unmapped))types = c("uniquely_mapped", "multimapped", "unmapped")to = ti2 %>% select(SampleID, uniquely_mapped, multimapped, unmapped)to
}

说明：

read_multiqc_bwa函数与read_multiqc_star函数相似，处理的是BWA（Burrows-Wheeler Aligner）比对的输出数据。BWA是另一个常用的比对工具，尤其适用于DNA-Seq数据。这个函数的主要任务是读取BWA的比对统计信息，并提取关键的比对数据。

读取数据：
- ti = read_tsv(fi)：读取MultiQC报告中的BWA比对统计数据。
处理成对和单端数据：
- 如果paired == F，则直接处理单端数据。
- 如果paired == T，则按样本ID分割，提取样本信息。
计算不同的比对情况：
- transmute(SampleID = SampleID, total = total_reads, ...)：计算总读取数、唯一比对数、多重比对数、未比对数等。
- nd = total - uniquely_mapped - multimapped - unmapped：计算丢失的读取数（nd）。
数据汇总：
- 使用group_by(SampleID)按样本ID分组，汇总每个样本的比对统计信息。
- 最终返回一个包含每个样本唯一比对数、多重比对数和未比对数的数据框。

7. 定义`read_multiqc_fastqc`函数

read_multiqc_fastqc <- function(fi) {ti = read_tsv(fi)ti2 = ti %>%mutate(SampleID = Sample) %>%select(-Sample)types = c("pass", "fail")to = ti2 %>% select(SampleID, total = total_sequences,passed = passed, failed = failed)to
}

说明：

read_multiqc_fastqc函数用于读取并处理FastQC质量控制工具生成的报告。FastQC是一个非常常用的工具，用于检查高通量测序数据的质量，如GC含量、序列长度分布等。此函数主要提取与数据质量相关的统计信息。

读取数据：
- ti = read_tsv(fi)：读取MultiQC报告中的FastQC数据。
处理数据：
- 直接将样本ID从Sample列提取，并将其命名为SampleID。
输出：
- 返回一个包含样本ID、总序列数（total_sequences）、通过质量控制的序列数（passed）和未通过质量控制的序列数（failed）的数据框。

8. 定义`get_input_reads_multiqc`函数

get_input_reads_multiqc <- function(multiqc_fi) {ti = read_multiqc_trimmomatic(multiqc_fi)to = ti %>% transmute(SampleID = SampleID, total = total,surviving = surviving) to
}

说明：

get_input_reads_multiqc函数用于从MultiQC报告中提取输入的读取数。它通过调用之前定义的read_multiqc_trimmomatic函数，提取Trimmomatic处理后的统计信息，并返回相关的数据框。

调用read_multiqc_trimmomatic：
- ti = read_multiqc_trimmomatic(multiqc_fi)：调用之前定义的read_multiqc_trimmomatic函数来读取并处理MultiQC报告中的Trimmomatic统计数据。
输出：
- 返回一个包含样本ID、总读取数和存活读取数的数据框。

9. 定义`get_align_summary_multiqc`函数

get_align_summary_multiqc <- function(multiqc_fi, aligner = "bwa") {if (aligner == "bwa") {ti = read_multiqc_bwa(multiqc_fi)} else if (aligner == "star") {ti = read_multiqc_star(multiqc_fi)} else {stop("unsupported aligner")}ti
}