拼接两条有重叠区域的核酸序列
目的
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 包完成：- 导入序列；
- 反补序列；
- 联配；
- 生成合并序列；
- 导出序列。
安装
步骤
1 准备
- 新建文件夹，放入：
- 代码
merge2Seqs.v1.R fasta格式的正反向序列。- seq1.fas
- seq2.fas
- 代码
不是要求文件后缀为
fasta，而是文件的内容符合每条序列第一行为>开头。
2 运行
- 打开
CMD。 - 进入文件夹所在路径。
cd your_path- 运行代码帮助
RScript merge2Seqs.v1.R -help- 运行代码
RScript merge2Seqs.v1.R seq1.fas seq2.fas titletitle如果不提供，默认为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 20181119
options(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 strand
er1 <- 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) <- title
fileout <- 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 = " "))
----ฅ(*ΦωΦ)ฅ---- cognata ad sidera tendit...
