

R自定义函数进行核酸序列格式转换

R语言数据分析指南

2023-03-25

欢迎关注R语言数据分析指南

❝
之前在处理生信数据的时候有时会有将核酸序列转换成氨基酸序列的需求，之前总是借助各种已有的工具来进行数据的转换，近来由于ChatGPT越来越强大于是借助其来实现这样功能。下面来进行代码介绍过程，「代码整理已经上传2023VIP群，加群的观众老爷有兴趣可请自行下载运行」
❞

❝
有需要学习个性化数据可视化的朋友，欢迎到小编的「淘宝店铺」 「R语言数据分析指南」购买「2023年度会员」 初始价格99元，内容主要包括各种「高分论文的图表分析复现以及一些个性化图表的绘制」均包含数据+代码；按照往年数据小编年产出约在150篇以上。
❞

购买后微信发小编订单截图即邀请进新的会员交流群，小编的文档为按年售卖，只包含当年度的「除系列课程外」的文档，有需要往年文档的朋友也可在淘宝店铺下单购买，需要了解更多信息的朋友欢迎交流咨询，「添加小编微信请备注，以提高沟通效率」

加载R包

library(tidyverse)

建立核苷酸和氨基酸的字典

nt_dict <- list(TTT="F", TTC="F", TTA="L", TTG="L",
                CTT="L", CTC="L", CTA="L", CTG="L",
                ATT="I", ATC="I", ATA="I", ATG="M",
                GTT="V", GTC="V", GTA="V", GTG="V",
                TCT="S", TCC="S", TCA="S", TCG="S",
                CCT="P", CCC="P", CCA="P", CCG="P",
                ACT="T", ACC="T", ACA="T", ACG="T",
                GCT="A", GCC="A", GCA="A", GCG="A",
                TAT="Y", TAC="Y", TAA="*", TAG="*",
                CAT="H", CAC="H", CAA="Q", CAG="Q",
                AAT="N", AAC="N", AAA="K", AAG="K",
                GAT="D", GAC="D", GAA="E", GAG="E",
                TGT="C", TGC="C", TGA="*", TGG="W",
                CGT="R", CGC="R", CGA="R", CGG="R",
                AGT="S", AGC="S", AGA="R", AGG="R",
                GGT="G", GGC="G", GGA="G", GGG="G")

将核苷酸序列转换为氨基酸序列

translate_dna <- function(dna_seq) {
  nts <- strsplit(dna_seq, "")[[1]]
  codons <- paste0(nts[c(T, F, F)], nts[c(F, T, F)], nts[c(F, F, T)])
  aa_seq <- sapply(codons, function(codon) nt_dict[[codon]])
  return(paste(aa_seq, collapse=""))
}

测试代码

dna_seq <- "ATGGTCTGCTGTGAC"
translate_dna(dna_seq)

❝
上方代码可以将核酸序列转换成氨基酸序列，但是如果我们拥有多段序列那就需要批量进行处理，下面介绍两种方式批量进行序列的翻译。
❞

Biostrings进行序列翻译

BiocManager::install("Biostrings")
library(Biostrings)
seqs <- readDNAStringSet("HvOSCA.cds.fa")

翻译核酸序列为氨基酸序列

aa_seqs <- translate(seqs)
writeXStringSet(aa_seqs, "output_aa_seqs.fasta")

❝
可以看到使用R包非常的方便，极大的简化了代码，下面让我们来尝试自定义代码来进行序列翻译
❞

自定义函数进行序列翻译

translate_fasta_dna_to_aa <- function(file_in, file_out) {
  # 读取FASTA文件
  fasta_file <- readLines(file_in)
  
  # 从FASTA文件中提取序列
  seqs <- list()
  for (i in 1:length(fasta_file)) {
    if (startsWith(fasta_file[i], ">")) {
      header <- gsub(">", "", fasta_file[i])
      seqs[[header]] <- ""
    } else {
      seqs[[header]] <- paste0(seqs[[header]], fasta_file[i])
    }
  }
  
  # 核酸翻译表
  codon_table <- list(
    "TTT" = "F", "TTC" = "F", "TTA" = "L", "TTG" = "L",
    "TCT" = "S", "TCC" = "S", "TCA" = "S", "TCG" = "S",
    "TAT" = "Y", "TAC" = "Y", "TAA" = "*", "TAG" = "*",
    "TGT" = "C", "TGC" = "C", "TGA" = "*", "TGG" = "W",
    "CTT" = "L", "CTC" = "L", "CTA" = "L", "CTG" = "L",
    "CCT" = "P", "CCC" = "P", "CCA" = "P", "CCG" = "P",
    "CAT" = "H", "CAC" = "H", "CAA" = "Q", "CAG" = "Q",
    "CGT" = "R", "CGC" = "R", "CGA" = "R", "CGG" = "R",
    "ATT" = "I", "ATC" = "I", "ATA" = "I", "ATG" = "M",
    "ACT" = "T", "ACC" = "T", "ACA" = "T", "ACG" = "T",
    "AAT" = "N", "AAC" = "N", "AAA" = "K", "AAG" = "K",
    "AGT" = "S", "AGC" = "S", "AGA" = "R", "AGG" = "R",
    "GTT" = "V", "GTC" = "V", "GTA" = "V", "GTG" = "V",
    "GCT" = "A", "GCC" = "A", "GCA" = "A", "GCG" = "A",
    "GAT" = "D", "GAC" = "D", "GAA" = "E", "GAG" = "E",
    "GGT" = "G", "GGC" = "G", "GGA" = "G", "GGG" = "G"
  )
  
  # 翻译核酸序列为氨基酸序列
  translate_seq <- function(dna_seq) {
    n_codons <- floor(nchar(dna_seq) / 3)
    codons <- substring(dna_seq, seq(1, n_codons * 3, 3), seq(3, n_codons * 3, 3))
    aa_seq <- ""
    for (codon in codons) {
      aa <- codon_table[[codon]]
      aa_seq <- paste0(aa_seq, aa)
    }
    return(aa_seq)
  }  
  #对每条序列进行翻译
  aa_seqs <- list()
  for (header in names(seqs)) {
    dna_seq <- seqs[[header]]
    aa_seq <- translate_seq(dna_seq)
    aa_seqs[[header]] <- aa_seq
  }
  
  # 将氨基酸序列输出到FASTA文件
  output_file <- file(file_out, "w")
  for (header in names(aa_seqs)) {
    aa_seq <- aa_seqs[[header]]
    cat(">", header, "\n", aa_seq, "\n", file = output_file, sep = "")
  }
  close(output_file)
}

使用案例

translate_fasta_dna_to_aa("HvOSCA.cds.fa", "output_aa_seqs.fasta")

❝
可以看到自定义函数进行处理还是比较麻烦的不如直接调用R包来的简便，下面再来介绍一下如何将氨基酸序列转换为核酸序列。
❞

氨基酸序列转核酸

library(Biostrings)
aa_seqs <- readAAStringSet("input.fasta") # 读入氨基酸序列文件
codon_table <- DNACodonTable(1) # 定义遗传密码表
# 将氨基酸序列翻译为核酸序列
dna_seqs <- translate(aa_seqs, codon_table)
# 将核酸序列输出到FASTA文件
writeXStringSet(dna_seqs, "output.fasta")

将一段序列进行翻译

# 定义一个核酸序列
dna_seq <- DNAString("ATGGGTTCTCTCAACGAAATCGGCGTTGCTGCGGGAATAAACATATCATCGGCATTGGGTTTTCTTCTAG")
#将核酸序列转换为AAString对象
aa_seq <- translate(dna_seq)
#将AAString对象转换为字符串
aa_seq_str <- as.character(aa_seq)
#输出氨基酸序列
cat(aa_seq_str)