R语言数据处理那些事（九） 蛋白序列覆盖度计算

代码

library(data.table)
library(tidyverse)
 
#考虑 测序的序列可能重叠，因此需要定义合并重叠肽段的函数
merge_segments <- function(segments) {
    # 如果只有一个肽段，直接返回
    if (nrow(segments) == 1) {
        return(segments)
    }
    
    merged <- list()
    current_segment <- segments[1, ]
    
    for (i in 2:nrow(segments)) {
        next_segment <- segments[i, ]
        
        # 如果当前肽段和下一个肽段有重叠或相邻
        if (current_segment$end >= next_segment$start) {
            # 合并这两个肽段
            current_segment <- data.frame(start = current_segment$start,
                                          end = max(current_segment$end, next_segment$end))
        } else {
            # 如果没有重叠，保存当前肽段，并开始新的合并
            merged <- append(merged, list(current_segment))
            current_segment <- next_segment
        }
    }
    
    merged <- append(merged, list(current_segment))
    
    # 返回合并后的结果
    return(do.call(rbind, merged))
}
 
#定义计算覆盖度的函数
cal_coverage <- function(peptides, sequence){
    # 合并重叠肽段
    pos <- lapply(peptides, function(x) str_locate_all(sequence, x)[[1]])
    # 转换为data.table，便于快速处理
    pos_dt <- rbindlist(lapply(pos, function(p) data.table(start = p[, 1], end = p[, 2])))
    
    # 排序并合并重叠肽段
    pos_sorted <- pos_dt[order(start)] # 排序
    merged_pos <- merge_segments(pos_sorted)
    
    # 计算覆盖度
    coverage <- round((sum(merged_pos$end - merged_pos$start) + nrow(merged_pos)) / nchar(sequence), 2)
    
    return(coverage)
}
 
#测试数据源于maxquant软件结果测试函数是否正确
sequence1 <- "MATVAATTKVPEIRDVTRIERIGAHSHIRGLGLDDALEPRQASQGMVGQLAARRAAGVVLEMIREGKIAGRAVLIAGQPGTGKTAIAMGMAQALGPDTPFTAIAGSEIFSLEMSKTEALTQAFRRSIGVRIKEETEIIEGEVVEIQIDRPATGTGSKVGKLTLKTTEMETIYDLGTKMIESLTKDKVQAGDVITIDKATGKISKLGRSFTRARDYDAMGSQTKFVQCPDGELQKRKEVVHTVSLHEIDVINSRTQGFLALFSGDTGEIKSEVREQINAKVAEWREEGKAEIIPGVLFIDEVHMLDIESFSFLNRALESDMAPVLIMATNRGITRIRGTSYQSPHGIPIDLLDRLLIVSTSPYSEKDTKQILRIRCEEEDVEMSEDAYTVLTRIGLETSLRYAIQLITAASLVCRKRKGTEVQVDDIKRVYSLFLDESRSTQYMKEYQDAFLFNELKGETMDTS"
peps1 <- c("AAGVVLEMIR", "AVLIAGQPGTGK", "EVVHTVSLHEIDVINSR","EYQDAFLFNELK","GLGLDDALEPR", "GTSYQSPHGIPIDLLDR","VYSLFLDESR") 
 
 
sequence2 <- "MREIVHLQAGQCGNQIGAKFWEVISDEHGIDPTGTYHGDSDLQLERINVYYNEATGGKYVPRAVLVDLEPGTMDSVRSGPFGQIFRPDNFVFGQSGAGNNWAKGHYTEGAELVDSVLDVVRKEAESCDCLQGFQLTHSLGGGTGSGMGTLLISKIREEYPDRIMNTFSVVPSPKVSDTVVEPYNATLSVHQLVENTDETYCIDNEALYDICFRTLKLTTPTYGDLNHLVSATMSGVTTCLRFPGQLNADLRKLAVNMVPFPRLHFFMPGFAPLTSRGSQQYRALTVPELTQQMFDAKNMMAACDPRHGRYLTVAAVFRGRMSMKEVDEQMLNVQNKNSSYFVEWIPNNVKTAVCDIPPRGLKMSATFIGNSTAIQELFKRISEQFTAMFRRKAFLHWYTGEGMDEMEFTEAESNMNDLVSEYQQYQDATAEEEGEFEEEAEEEVA"
peps2 <- c("INVYYNEATGGK","EIVHLQAGQCGNQIGAK","FWEVISDEHGIDPTGTYHGDSDLQLER","ALTVPELTQQMFDAK","EAESCDCLQGFQLTHSLGGGTGSGMGTLLISK","IMNTFSVVPSPK","YLTVAAVFR","IREEYPDR","FPGQLNADLR","LAVNMVPFPR","EVDEQMLNVQNK","GHYTEGAELVDSVLDVVR","ISEQFTAMFR","LHFFMPGFAPLTSR","LTTPTYGDLNHLVSATMSGVTTCLR","SGPFGQIFRPDNFVFGQSGAGNNWAK","NMMAACDPR","NSSYFVEWIPNNVK","MSATFIGNSTAIQELFK") 
 
cal_coverage(peps1, sequence1) #0.19 和maxquant计算结果一致
cal_coverage(peps2, sequence2) #0.67 和maxquant计算结果一致