实验课3 第七周实验课 R代码整理¶

对应材料：

实验/课件/实验课3-第七周实验课_修改版.pptx
老师板书整理

说明¶

下面这份代码是按老师板书主线整理的“可直接复习版”。
板书里有几处只写了函数骨架、没有完全写完，我已按课件中的明确内容补齐。
补齐时尽量保持老师板书风格，不额外引入复杂写法。

一、板书整理后的完整 R 代码¶

library(Biostrings)

file_path <- "DNA_seq.fasta"
lines <- readLines(file_path, warn = FALSE)
seq_lines <- lines[!grepl("^>", lines)]
dna_seq <- toupper(paste(seq_lines, collapse = ""))

cat(">>> 序列长度:", nchar(dna_seq), "\n\n")

expand_inr_motifs <- function() {
  pos1 <- c("C", "T")
  pos2 <- c("C", "T")
  pos3 <- c("A")
  pos4 <- c("A", "T", "C", "G")
  pos5 <- c("A", "T")
  pos6 <- c("C", "T")
  pos7 <- c("C", "T")

  combs <- expand.grid(
    pos1, pos2, pos3, pos4, pos5, pos6, pos7,
    stringsAsFactors = FALSE
  )

  apply(combs, 1, paste0, collapse = "")
}

predict_TSS <- function(seq) {
  s <- DNAString(seq)
  tss_candidates <- data.frame()

  motifs <- list(
    TATA = c("TATAAAA", "TATAAAT", "TATATAA", "TATATAT"),
    Inr = expand_inr_motifs(),
    BRE = c("GGAGCC", "GGTGCC", "GCAGCC", "GCTGCC"),
    CAAT = c("GGCCAATCT"),
    GC = c("GGGCGG")
  )

  for (name in names(motifs)) {
    for (mot in motifs[[name]]) {
      matches <- matchPattern(DNAString(mot), s)

      if (length(matches) > 0) {
        tss_candidates <- rbind(
          tss_candidates,
          data.frame(
            Motif = name,
            Motif_seq = mot,
            Start = start(matches),
            End = end(matches),
            stringsAsFactors = FALSE
          )
        )
      }
    }
  }

  if (nrow(tss_candidates) == 0) {
    return(tss_candidates)
  }

  tss_candidates$Predicted_TSS <- NA_integer_

  idx <- tss_candidates$Motif == "TATA"
  tss_candidates$Predicted_TSS[idx] <- tss_candidates$End[idx] + 25

  idx <- tss_candidates$Motif == "BRE"
  tss_candidates$Predicted_TSS[idx] <- tss_candidates$End[idx] + 35

  idx <- tss_candidates$Motif == "CAAT"
  tss_candidates$Predicted_TSS[idx] <- tss_candidates$End[idx] + 100

  idx <- tss_candidates$Motif == "GC"
  tss_candidates$Predicted_TSS[idx] <- tss_candidates$End[idx] + 80

  idx <- tss_candidates$Motif == "Inr"
  tss_candidates$Predicted_TSS[idx] <- floor(
    (tss_candidates$Start[idx] + tss_candidates$End[idx]) / 2
  )

  tss_candidates
}

predict_TSE <- function(seq) {
  s <- DNAString(seq)
  tse_candidates <- data.frame()

  polyA_motifs <- c("AATAAA", "ATTAAA", "AGTAAA", "TATAAA")

  for (mot in polyA_motifs) {
    matches <- matchPattern(DNAString(mot), s)

    if (length(matches) > 0) {
      tse_candidates <- rbind(
        tse_candidates,
        data.frame(
          Motif = mot,
          Start = start(matches),
          End = end(matches),
          stringsAsFactors = FALSE
        )
      )
    }
  }

  if (nrow(tse_candidates) == 0) {
    return(tse_candidates)
  }

  tse_candidates$Predicted_TSE <- tse_candidates$End + 20
  tse_candidates
}

true_tss <- 2018
true_tse <- 194629

tss_candidates <- predict_TSS(dna_seq)
tse_candidates <- predict_TSE(dna_seq)

if (nrow(tss_candidates) > 0) {
  tss_candidates$TSS_error <- abs(tss_candidates$Predicted_TSS - true_tss)
  tss_candidates <- tss_candidates[order(tss_candidates$TSS_error), ]
}

if (nrow(tse_candidates) > 0) {
  tse_candidates$TSE_error <- abs(tse_candidates$Predicted_TSE - true_tse)
  tse_candidates <- tse_candidates[order(tse_candidates$TSE_error), ]
}

cat(">>> TSS candidates:\n")
print(head(tss_candidates, 10))

cat("\n>>> TSE candidates:\n")
print(head(tse_candidates, 10))

二、这份代码对应老师板书的哪几块¶

1. 数据读入与预处理¶

对应老师板书右上角：

library(Biostrings)
file_path <- "DNA_seq.fasta"
lines <- readLines(file_path, warn = FALSE)
seq_lines <- lines[!grepl("^>", lines)]
dna_seq <- toupper(paste(seq_lines, collapse = ""))

意思就是：

读取 FASTA
去掉 > 开头的标题行
拼接多行 DNA 序列
统一转成大写

2. `expand_inr_motifs()`¶

对应老师板书中间和左侧那块：

pos1 <- c("C", "T")
pos2 <- c("C", "T")
pos3 <- c("A")
pos4 <- c("A", "T", "C", "G")
pos5 <- c("A", "T")
pos6 <- c("C", "T")
pos7 <- c("C", "T")

再用：

expand.grid(...)
apply(..., paste0, collapse = "")

把所有可能的 Inr 序列展开出来。

这其实就是把课件里的：

YYAN[AT]YY

具体展开成所有可能组合。

3. `predict_TSS()`¶

对应第二张板书图：

先建立 motifs 列表
再双重循环：
外层循环：不同 motif 类型
内层循环：该 motif 类型下的每一条具体序列
用 matchPattern() 扫描序列
如果扫到匹配结果，就把结果存进 data.frame

板书里已经明确写出了这几个核心对象：

s <- DNAString(seq)
tss_candidates <- data.frame()
matches <- matchPattern(DNAString(mot), s)
start(matches)
end(matches)

4. `predict_TSE()`¶

对应第三张板书图：

先定义 polyA_motifs
逐个扫描
把命中的 Start、End 保存下来
再按老师板书直接写：

tse_candidates$Predicted_TSE <- tse_candidates$End + 20

三、我按课件补全了哪些地方¶

1. `TSS` 的具体 motif 序列¶

板书只写了 TATA、Inr、BRE、CAAT、GC 这几个名字，没有把每一条序列都写全。

这里我按课件补成：

TATA = c("TATAAAA", "TATAAAT", "TATATAA", "TATATAT")
BRE  = c("GGAGCC", "GGTGCC", "GCAGCC", "GCTGCC")
CAAT = c("GGCCAATCT")
GC   = c("GGGCGG")

2. `TSS` 的预测位置计算¶

板书第二张图主要写到“候选位点收集”这一步，还没完全把 Predicted_TSS 那几行展开。

这里我按课件规则补成：

TATA：下游 +25
BRE：下游 +35
CAAT：下游 +100
GC：下游 +80
Inr：取中心位置

为了和老师板书中 TSE = End + 20 的写法保持一致，这里我统一写成了“以匹配片段末端 End 为基准，再加偏移量”。

3. 最后的误差分析¶

板书右下角写了：

TSS <- 2018
TTS <- 194629

这里课件和我们前面整理统一写作：

true_tss <- 2018
true_tse <- 194629

然后补上误差分析：

abs(Predicted_TSS - true_tss)
abs(Predicted_TSE - true_tse)

这样就能直接看哪个候选位点更接近真实值。

四、你复习时最该盯住的代码逻辑¶

如果考试不要求你逐行默写代码，那最值得记的是这条主线：

读入FASTA
-> 去标题行并拼接序列
-> 扫描TSS相关motif
-> 保存所有候选TSS
-> 扫描polyA motif
-> 保存所有候选TSE
-> 和真实位点比较误差

其中最容易单独拿出来考的，是这四个点：

FASTA 怎么预处理
Inr 为什么要展开组合
双重循环在什么地方
data.frame 为什么要保存全部候选位点

五、如果你考试要手写伪代码，可以写成这样¶

输入：一条FASTA格式DNA序列

1. 读取FASTA文件
2. 去掉以 > 开头的标题行
3. 拼接序列并转为大写
4. 定义TSS相关motif：TATA、Inr、BRE、CAAT、GC
5. 对每一类motif逐个扫描整条序列
6. 如果发现匹配，则记录motif类型、起点、终点
7. 按不同motif规则换算出候选TSS
8. 定义polyA motifs
9. 再次扫描整条序列
10. 如果发现匹配，则记录起点、终点，并换算出候选TSE
11. 将候选结果与真实TSS/TSE比较误差
12. 输出候选位点和误差结果

六、一句提醒¶

这份代码最重要的不是“函数名必须一字不差”，而是你要会解释：

输入是什么
循环在哪里
判断在哪里
输出是什么
为什么这样能从 DNA 序列里初步预测 TSS/TSE