基因组Survey分析

流程图：
在这里插入图片描述
图片来源：https://www.jianshu.com/p/94da86093843

一、Fastp质控

在这里插入图片描述

二、NT比对

一般选择第六个输出格式
在这里插入图片描述

结果示例：
在这里插入图片描述

三、k-mer分析

软件：GCE/genomescope
分析目的：预估基因组大小，重复序列比例
k-mer：核酸序列分成包含K个碱基的字符串
mer代表monomeric unit，可看成单bp。
核酸序列长度为L，k-mer长度为K，那么可以得到L-K+1个k-mer

在这里插入图片描述
使用k-mer估计基因组大小：

一般以主峰深度作为kmer期望测序深度，也就是那个分号下面的
实操：
获取测序数据的K-mer频率：

kmer -k 17 -t 20 -p qianzhui cleandatalist

-k 指设置的K值，一般为17 一定要在13-19之间，不然会报错
-p 指输出文件前缀
cleandatalist 是一个包含路径的文件，如下：
在这里插入图片描述
生成名为kmer.freq.stat的文件，前缀是qianzhui，所以文件就是qianzhui.kmer.freq.stat，使用这个文件进行下一步分析
本步骤参考：https://blog.csdn.net/weixin_69890544/article/details/135440067

获取gce运行参数：

less qianzhui.kmer.freq.stat | grep "#Kmer indivdual number" #获取gce参数-g
less qianzhui.kmer.freq.stat | perl -ne 'next if(/^#/ || /^\s/); print; ' | awk '{print $1"\t"$2}' > ara.kmer.freq.stat.2colum 
#获得gce参数-f，也就是这里的ara.kmer.freq.stat.2colum

使用gce进行survey：
纯合模式：

gce -g 3295248520 -f ara.kmer.freq.stat.2colum >gce.table 2>gce.log
#使用之前的得到的-g和-f参数进行基因组survey

杂合模式：

./gce -g 3295248520 -f ara.kmer.freq.stat.2colum -H 1 -c 28 >gce.table 2>gce.log
#使用之前的得到的-g和-f参数进行基因组survey
# -c 期望深度 其中-c 的值一般指定为纯合模式运行得到的rawpeak，在纯合得到的gce.log里面有

生成的文件里有一行是kmer-species heterozyugous ratio is
后面有一个数字，用这个数字除以kmer大小，等于基因组杂合率，基因组杂合率<0.002,可以大概判断是纯合基因组，否则是杂合基因组
根据是纯合还是杂合基因组，去使用对应的gce.log的结果
genomesize，是基因组大小。
重复序列占比：
纯合模式运行结果的最下面有genomesize和b[1]信息
杂合模式的运行结果最下面有genomesize和b[1],b[1/2]信息
纯合模式重复序列占比=1-b[1]
杂合模式重读序列占比=1-b[1/2]-b[1]
本步骤参考：https://blog.csdn.net/weixin_69890544/article/details/135440067

四、画个图

c=`awk '$1==60' ara.kmer.freq.stat.2colum|awk '{print $2}'`
echo $c
#选取合理的深度范围
head -n 500 ara.kmer.freq.stat.2colum > ara.freq.stat.2colum.500
#作图
Rscript distribution.r ara.freq.stat.2colum.500 ./ $c
convert kmer_distribution.svg kmer_distribution.png
sz kmer_distribution.png

distribution.r 脚本内容：

library(ggplot2)
#1. data # 读入 深度-Kmer种类数频率 表格
args <- commandArgs()
file=args[6] 
a<-read.table(file,sep="\t")
#2. output
setwd(args[7])
#3. ylim 峰值大小，就是Kmer的种类数峰值大小，作为y的max值
peak=args[8]
peak<-as.numeric(peak)
#4. plot 作图，
svg("kmer_distribution.svg", width=10) 
ggplot(a,aes(x=V1,y=V2),col="red")+geom_line(color="green")+geom_point(color="red")+xlim(0,200)+ylim(0,peak)+xlab("
Depth of Kmer Species")+ylab("Frequency of Kmer Species")+theme_bw()+theme(axis.title=element_text(size=20))
dev.off()