ls /Data/A* | awk '{ print "gzip -dc " $0}' > A_generate.file
jellyfish count -t 20 -C -m 21 -s 3G -g A_generate.file -G 2 -o A.sif
# jellyfish count:k-mer计数
# -t 指定线程数
# -C 对DNA正负链都进行统计,表示考虑DNA正义与反义链。如果是双端测序reads,需要这个参数。
# -m k-mer长度设置为21bp.
# -s 3G 存储用的hash表大小为3G,这个参数识别单位M(Mbp)和G(Gbp)。如果设置不够大,会生成多个hash文件。最好设置的值大于总的独特的(distinct)k-mer数。如果基因组大小为G,每个reads有一个错误,总共有n条reads,则该值可以设置为[(G + n)/0.8]。
# -g --generator=path 记录产生fast[aq]命令的文件
# -G 同时运行的数目
# -o 指定输出文件的名字,hash格式储存的k-mer频数文件
jellyfish histo -v -t 20 -h 10000000 -o A.histo A.sif
# 结果最后是两列,第一列是x,表示的是出现的次数;第二列是y,表示出现x次数的kmer数目。
# -v 显示详细信息
# -t 线程数
# -h x的最大值,默认是10000
# -o 输出文件的名字。
# 估计覆盖阈值上(U)下(L)限
U=$(smudgeplot.py cutoff A.histo U)
L=$(smudgeplot.py cutoff A.histo L)
# 计算k-mer pairs
jellyfish dump -c -L $L -U $U A.sif | smudgeplot.py hetkmers -o kmer_pairs
# dump:将Jellyfish计数后的k-mer数据转储(dump)出来。
# -c 参数表示输出覆盖度(coverage)信息
# -L $L 和 -U $U 参数分别表示覆盖度的下限(lower bound)和上限(upper bound),确定的合理阈值。
# smudgeplot.py hetkmers 寻找杂合k-mer pairs
# 绘图
smudgeplot.py plot kmer_pairs_coverages.tsv -o my_genome
# 结果图片:my_genome_smudgeplot.png
# 具体数值:my_genome_summary_table.tsv
mkdir tmp
ls /Data/A* | awk '{ print "gzip -dc " $0}' > A_generate.file
kmc -k21 -t16 -m64 -ci1 -cs10000 A_generate.file kmcdb tmp #计算k-mer频率
# -k21:k-mer长度设置为21
# -t 指定线程数
# -m64:内存64G,设置使用RAM的大致数量,范围1-1024。
# -ci1 -cs10000:统计k-mer coverages覆盖度范围在[1-10000]。
# A_generate.file:保存了输入文件列表的文件名为A_generate.file
# kmcdb:KMC数据库的输出文件名前缀
# tmp:临时目录
kmc_tools transform kmcdb histogram A.histo -cx10000 #生成k-mer频数直方表A.histo和k-mer直方图
# -cx10000:储存在直方图文件中counter的最大值。
# 估计覆盖阈值上(U)下(L)限
U=$(smudgeplot.py cutoff A.histo U)
L=$(smudgeplot.py cutoff A.histo L)
# 用kmc_tools提取k-mers
kmc_tools transform kmcdb -ci"$L" -cx"$U" reduce kmcdb_L"$L"_U"$U"
# 用KMC的smudge_pairs计算k-mer pairs
smudge_pairs kmcdb_L"$L"_U"$U" kmcdb_L"$L"_U"$U"_coverages.tsv kmcdb_L"$L"_U"$U"_pairs.tsv > kmcdb_L"$L"_U"$U"_familysizes.tsv
# 绘图
smudgeplot.py plot kmcdb_L"$L"_U"$U"_coverages.tsv -o my_genome
# 结果图片:my_genome_smudgeplot.png
# 具体数值:my_genome_summary_table.tsv
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!