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)