目的

• 27F 和 1492R 是细菌 16S 核糖体基因的一对通用引物。利用这对引物扩增DNA提取物，然后进行Sanger法测序，能获得两条长度约为 800+ 的 DNA 序列。根据正链（+）和负链（-）序列之间的重叠部分拼接，就能拼接出一条 16S rRNA 基因序列。
  • F = forward，R = Reverse，数字 = 结合的位置。
  • 27F： TACGGYTACCTTGTTACGACTT。
  • 1492R: AGAGTTTGATCMTGGCTCAG。
• 因为序列书写和合成都是认从 5-P端到3'-OH端 ，所以（-）链反向之后，才能与（+）链互补。而要从（-）推出对应的（+）链部分，需要进行 反补（Reverse complement）。
• 拼接 的原理：利用 局部联配（local alignment） 找到 重叠的区域（overlap）。

从测序到拼接过程示意图

引物与序列

• DNA双链模板：

  5' ATGCAGGGGAAACATGATTCAGGAC 3' （+）
  3' TACGTCCCCTTTGTACTAAGTCCTG 5' （-）

• 正向引物 Forward primer：ATG；

• 反向引物 Reverse primer：GTC。

引物与序列结合

• F引物 ATG，与（-）链结合，ATG 符合5'到3'的方向，延伸得到 （+）链。
• R引物 GTC，与（+）链结合，G在5'端，从G往C合成 ，得到 （-）链。

                           CTG(R)
  5' ATGCAGGGGAAACATGATTCAGGAC 3' （+）
  3' TACGTCCCCTTTGTACTAAGTCCTG 5' （-）ATG(F)

测序结果

• 按惯例，序列写作从5'端到3'端。

• （+）链还是正的：

  5' ATGCAGGGGAAACATGATTCAGGAC 3' （+）模板
  5' ATGCAGGGGAAACATG————————— 3' （+）合成结果

• （-）链的结果反向书写：

  3' TACGTCCCCTTTGTACTAAGTCCTG 5' （-）模板
  5' GTCCTGAATCATGTTTCCCCTGCAT 3' （-）模板反向
  5' GTCCTGAATCATGTTT————————— 3' （-）合成结果

• （-）链反补推出（+）链部分：

  5' GTCCTGAATCATGTTT————————— 3' （-）合成结果
  3' —————————TTTGTACTAAGTCCTG 5' （-）反向
  5' —————————AAACATGATTCAGGAC 3' （-）反向，碱基互补 = （+）

• 根据重叠区域拼接，也相当于局部联配，Overlap = local alignment。

  5' ATGCAGGGGAAACATG————————— 3' （+）合成结果
  5' —————————AAACATGATTCAGGAC 3' （-）反向，碱基互补 = （+）
  5' ATGCAGGGGAAACATGATTCAGGAC 3' （+）合并结果 consensus

工具

软件

• 使用 R 中的 DECIPHER 包完成：
  • 导入序列；
  • 反补序列；
  • 联配；
  • 生成合并序列；
  • 导出序列。

安装

• 安装 R > 3.5.0;
• RScript 能在 CMD 中运行 = R的环境变量设置正确。
• R 中已经安装 DECIPHER。
• 详细安装步骤请参考 link

步骤

1 准备

• 新建文件夹，放入：
  • 代码 merge2Seqs.v1.R
  • fasta 格式的正反向序列。
    • seq1.fas

    • seq2.fas

      不是要求文件后缀为 fasta，而是文件的内容符合每条序列第一行为 > 开头。

2 运行

1. 打开 CMD 。
2. 进入文件夹所在路径。

   cd your_path

3. 运行代码帮助

   RScript merge2Seqs.v1.R -help

4. 运行代码

   RScript merge2Seqs.v1.R seq1.fas seq2.fas title

title如果不提供，默认为consensus sequence。

结果说明

• 文件夹中生成合并后的序列.
• 屏幕显示联配的详情，注意检查重叠区域的长度，错配的个数，序列总长度，结果中是不是有简并碱基。通常错配的位置在个位数。
• 只有两条序列无法确认哪一条是正确的，
  • 如果是通过向一条链插入空白来配对，则合并的序列中包括这个位点，

  • 如果两条序列上的位点无法匹配，该序列里面会用简并碱基表示。
    例如：

    Seq1: ATC-A ATCGA
    Seq2: ATCCA ATCCA
    Con: ATCCA ATCSA

  • 简并碱基表示方法：

  A : A
  C : C
  G : G
  T : T
  M : AC
  R : AG
  W : AT
  S : CG
  Y : CT
  K : GT
  V : ACG
  H : ACT
  D : AGT
  B : CGT
  N : ACGT

• 如果由于序列质量导致错配 mismatch 过多，需要自行调节原始序列的裁剪范围，重新拼接。
  • 通常Sanger测序法结果的前15-40bp的质量较差，全长约为700-900bp。
  • 但是根据实际经验，即使完全不掐头去尾，也能拼出错配低至1的序列，且与基因组中的16SrRNA基因序列完全一致。所以，一定检查配对情况！

更新

• 添加了代码输出鉴别碱基的数量。
• 修改了gapopen = -5 gapextension = -2 在 AlignSeqs() 函数部分
  • 发现如果不改会出现和 pairwiseAlignment 结果不一样。
  • 合并2W06 的正反向结果时发现，它有10个错配（pairwiseAlignment），结果里有10个简并碱基（AlignSeqs()），但一个是插入空白，一个是插入空白加错配。
  • 但是这个例子还是很奇怪，最后虽然中间是对的但是两头比 27F 和 1492R 延伸得还要多 1608，中间却又是正确的。
  • 使用两个因为联配后发现的确是引物前后的序列也被测出来了。没有原始数据，所以不知道这是怎么来的，基因组里面的全长16S也只有 1525bp。

下一步

• 因为使用了两个联配，一个用来生成合并序列，一个用来输出联配结果，而我现在的知识无法保证两个联配函数的行为一致，我最好能自己把 AlignSeqs() 的结果输出出来。

代码

• 新建文本文件（.txt）改名为 merge2Seqs.v1.R ，复制粘贴下面的代码。注意系统是否显示后缀。

# @author Zheng Weishuang
# @purpose 合并两条正反向测序的序列
# @date 20181119options(encoding = "UTF-8")# 传递参数
args = commandArgs(trailingOnly=TRUE)if (as.numeric(R.Version()$major) < 3 | as.numeric(R.Version()$minor) <5) {cat("\n======================================\n(1)R version currently loaded is older than 3.5.0\nCheck the PATH, or upgrade R.======================================\n\n")
}# 无参数传递时提醒使用帮助
if(length(args) == 0){cat("\n======================================\nType Rscript merge2Seqs.R -help for help\n======================================\n\n")stop()
}# 运行帮助
if(args[1] == '-help'){cat("\n======================================\n(1)R version > 3.5.0\n(2)Open CMD (Windows) or Terminal (Mac)，loacate the path for merge2Seqs\n(3)Script eample: Rscript --vanilla merge2Seqs.R file1 file1 title\n(4)sequences must be in fasta format\n======================================\n\n")stop()
}# 如果在当前文件夹找不到序列
if(!file.exists(args[1])|!file.exists(args[2])){cat("\n======================================\nCan not find the sequences. \nWrong path or wrong file names。\n======================================\n\n")stop()
}# 加载 DECIPHER
file <- try(suppressPackageStartupMessages(library(DECIPHER)))# file <- try(library(DECIPHER))
if(inherits(file, "try-error")) {cat("\n======================================\nFail to load DECIPHER\n======================================\n\n")stop()}# function
merge2Seqs <- function(forwardPath = NA, reversePath = NA, title = NA){# the out come is the positive strander1 <- try(seq1 <- readDNAStringSet(forwardPath))if(inherits(er1, "try-error")) {cat("\n======================================\nFile1 is not in fasta format.\nThe first line should starts with '>'.\n")stop()} else {cat("\n======================================\nFile1 loaded\n")}er2 <- try(seq2 <- readDNAStringSet(reversePath))if(inherits(er2, "try-error")) {cat("\nfile2 is not in fasta format\nThe first line should starts with '>'.======================================\n")stop()} else {cat("\nfile2 loaded\n======================================\n\n")}seqs <- c(seq1, reverseComplement(seq2))
# 联配  align <- AlignSeqs(seqs, verbose = FALSE, gapOpening = -5,gapExtension = -2)# 合并的序列res <- ConsensusSequence(align, ignoreNonBases = T)# 序列名字names(res) <- "consensus sequence"fileout <- "consensus sequence.fas"if(nchar(title) > 0){names(res) <- titlefileout <- paste0(title,".fas")} # 导出序列writeXStringSet(res, fileout)# 验证结果  writePairwiseAlignments(pairwiseAlignment(seqs[[1]],seqs[[2]], type='local',gapOpening = 5, gapExtension = 2), stdout())cat("======================================\nCheck the overlap!\n")mismatch <- substring(as.character(res),1:nchar(as.character(res)),1:nchar(as.character(res))) # 把res里的结果从1到1，2到2，3到3，每个提取出来为一个if(length(unique(mismatch))!=4){cat("Watch out the ambiguity bases.\n");for(i in 1:length(IUPAC_CODE_MAP)[1]){cat(names(IUPAC_CODE_MAP)[i],":",IUPAC_CODE_MAP[i],'\n\n')}}print((mismatch[!grepl('[A|T|C|G]', names(mismatch))]))cat("\ngapOpening = -5,gapExtension = -2\nThe consensus sequence has been saved.\n======================================\n")return(res)
}# 运行拼接函数
merge2Seqs(forwardPath = args[1], reversePath = args[2], title = paste(args[c(-1,-2)], collapse = " "))