美文网首页基因组组装生信学习
基因组分析 K-mer 第0回 随机生成fasta文件

基因组分析 K-mer 第0回 随机生成fasta文件

作者: Jason数据分析生信教室 | 来源:发表于2021-06-04 10:17 被阅读0次

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)

相关文章

网友评论

    本文标题:基因组分析 K-mer 第0回 随机生成fasta文件

    本文链接:https://www.haomeiwen.com/subject/sbbxrltx.html