写在前面
- 以下内容均来自我在菲沙基因(Frasergen)暑期生信培训班上记录的课堂笔记
1.Kmer定义
Kmer定义- 杂合Kmer 杂合Kmer
- 杂合Kmer峰 杂合Kmer峰
- 重复Kmer峰 重复Kmer峰
- 杂合变化模拟Kmer图 杂合变化模拟Kmer图
2.基因组大小预估方式
预估基因组大小公式- 也就是基因组大小=(基因切割的Kmer数目)/(主峰深度)
- 杂合和重复会影响统计,而一般的基于Kmer预估基因组的软件会对此做处理,让结果更贴近真实。
3.Kmer分析实战
- 软件:GCE
3.1 下载安装gce软件
wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.2.tar.gz
tar -zxvf gce-1.0.2.tar.gz
3.2 使用gce软件中的kmerfreq脚本切割kmer并统计频率深度表格
gce-1.0.2/kmerfreq -k 17 -t 10 -p freq list_of_clean
#-k 17:切割kmer的长度
#list_of_clean是质控后的文件名
#cat list_of_clean
#/local_data1/pop_clean_1P.fastq.gz
#/local_data1/pop_clean_2P.fastq.gz
less freq.kmer.freq.stat|perl -ne 'next if(/^#/ || /^\s/); print; ' | awk '{print $1"\t"$2}' > freq.stat.2colum
#total kmer number, i.e. total number of kmer individuals
c=`less freq.kmer.freq.stat| grep "#Kmer indivdual number"|awk '{print $4}'`
echo $c
3.3 将上一行的$c也就是总的kmer种类数目和freq.stat.2colum也就是频率深度表格输入,得到基因组大小,杂合,重复信息 ,-H和-c的选择很重要
#纯合模式
gce -g $c -f freq.stat.2colum 2>gce.log
#杂合模式
gce -H 1 -g $c -c 60 -f freq.stat.2colum 2>heterozgyousgce.log
-g:kmer总数, 从kmerfreq分析结果获取
-f freq.stat.2colum:Kmer频率分布,从kmerfreq分析结果获取
-H 1:是否启动杂合模式(1是杂合模式,推算出杂合率, 0是非杂合模式没有杂合度)
-c 60: Kmer主峰深度,由gce自己选,或者根据情况自己选择(峰的选择很重要)
- GCE输出结果说明
3.4 Kmer作图
c=`awk '$1==60' freq.stat.2colum|awk '{print $2}'`
echo $c
#选取合理的深度范围
head -n 500 freq.stat.2colum > freq.stat.2colum.500
#作图
Rscript distribution.r freq.stat.2colum.500 ./ $c
convert kmer_distribution.svg kmer_distribution.png
sz kmer_distribution.png
- R作图脚本(上面的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()
- 结果
总结
- Survey分析内容回顾 Survey分析内容回顾
- Tips Tips
网友评论