0604 Rain
前言
关于k-mer的基础定义和用法，网上有太多的文章解释，什么估计基因长度，杂合程度，最优拼接...但很多时候道理都懂，就是过不好这一生。看了半天，还是回到了起点，所以k-mer到底是个什么鬼。在这里，我准备通过实践梳理一下k-mer的定义。当然不是纯理论，而是用R来虚拟一些序列，然后通过动手操作，来理解k-mer的内含。最终的目的是学会使用和理解去年(2020年)Nature Comunication出的两款基于K-mer的强力基因组分析软件Smudgeplot和Genome Scope2。
第0回 随机生成fasta文件
之所以是第0回，是因为这次还不会接触到k-mer的具体概念，而是为了之后的k-mer分析作准备而生成一些序列。
0.1 随机生成50bp的序列
各碱基的概率为A:20%, C: 30%, G:30%, T: 20%
library(Biostrings)
param_len_ref <- 50 #指定序列长度
narabi <- c("A","C","G","T") #指定序列种类
param_composition <- c(20, 30, 30, 20) #(A,C,G,T的顺序)指定各个碱基的出现概率
#开始生成序列
ACGTset <- rep(narabi, param_composition) #根据narabi中の4个碱基以及指定其出现的概率生成字符串
reference <- paste(sample(ACGTset, param_len_ref, replace=T), collapse="")
#从刚才的字符串ACGTset里进行可重复随机取样拼接起来，并把顺序打乱
reference #这样一段随机生成的序列就出现了
[1] "CATTTTTACAACGAGGATGACCTGTACTCGACGACCCGTGCCAGTAAAAG"
0.2 随机生成序列并保存成fasta格式
这次要生成70bp的序列，碱基构成为A: 23%, C:27%, G:28%, T:22%
out_f <- "hoge2.fasta" # 先给最后要输出的文件给起个名字
param_len_ref <- 70
narabi <- c("A","C","G","T")
param_composition <- c(23, 27, 28, 22)
param_desc <- "contig01" #FASTA文件的Description赋值
library(Biostrings)
- 生成序列同0.1
ACGTset <- rep(narabi, param_composition)
reference <- paste(sample(ACGTset, param_len_ref, replace=T), collapse="")
- 转换成FASTA格式
fasta <- DNAStringSet(reference)
names(fasta) <- param_desc #给FASTA添加description
fasta
DNAStringSet object of length 1:
width seq names
[1] 70 ATCCAGCCTCGGGCCGAACGTGAGCTGAG...AGTACTAAACATGGTCTTCTCCCCTCCC contig01
- 文件保存
writeXStringSet(fasta, file=out_f, format="fasta", width=50)
0.3 生成有多段序列的fasta
- 一开始的操作都一样
out_f <- "hoge3.fasta"
param_len_ref <- c(24, 103, 65, 49) #指定三条序列的长度
narabi <- c("A","C","G","T")
param_composition <- c(26, 27, 24, 23) # 同上
param_desc <- "contig"
- 生成序列
library(Biostrings)
ACGTset <- rep(narabi, param_composition)
hoge <- NULL
for(i in 1:length(param_len_ref)){ # 用loop给三条序列赋值
hoge <- c(hoge, paste(sample(ACGTset, param_len_ref[i], replace=T), collapse=""))
}
- 变换fasta格式
fasta <- DNAStringSet(hoge)
names(fasta) <- paste(param_desc, 1:length(hoge), sep="_")#添加description行的记述
fasta
- 文件保存
writeXStringSet(fasta, file=out_f, format="fasta", width=50)
