基因组survey

几种做基因组survey的方法
包括jellyfish、genomescope2.0、gce等软件
1.jellyfish计算k-mer频率

jellyfish count -C -m 21 -s 1000000000 -t 10 -o reads.jf *.fastq
注意：jellyfish只支持未压缩的fasta格式和fastq格式
######################################################
# -s 内存                                             
# -t 线程数                                           
# -m kmer值 如果覆盖率低或者错误率高，需要缩小kmer值     
# -C 统计正负链                    
######################################################
#计算k-mer频率的直方图 生成read.histo文件
$ jellyfish histo -t 10 reads.jf > reads.histo

结果可视化 R语言结果可视化

>histo <- read.table("readsSH.histo",header = F, sep = " ", stringsAsFactors = F)
>
>P=ggplot(histo, aes(x=V1, y=V2))+geom_line(size=0.8,color="#6495ed")+xlim(0,200)+
ylim(0,58000000)+theme_classic()+xlab("K-mer depth")+ylab("Kmer Freqence")+
theme(axis.title.x=element_text(size=15),axis.title.y=element_text(size=15),axis.text.y=element_text(size=12),axis.text.x=element_text(size=12))
>P
#在图中添加一条虚线
>P+geom_vline(aes(xintercept=51), colour="#BB0000", linetype="dashed",size=0.8)
#xintercept=51 该值为峰值的位置

Kmer depth-Freqence.jpg

2.genomescope2.0
软件安装

git clone https://github.com/tbenavi1/genomescope2.0.git
cd genomescope2.0/
mkdir ~/R_libs
echo "R_LIBS=~/R_libs/" >> ~/.Renviron
Rscript install.R

基因组大小和杂合度评估

genomescope.R

genomescope.R -i read.histo -o OUTPUT -k 21 -n name -p 倍性 
-i 输入jellyfish产生的read.histo文件
-o 输出文件夹
-k kmer值
-p 倍性 填1、2、3、4、5、6
-n 名字
-l 平均k-mer覆盖度的初始猜测
-m the cutoff for excluding high frequence k-mers from the analysis.分析中最高的k-mer分布频率

结果

GenomeScope
网页版本：http://qb.cshl.edu/genomescope/genomescope2.0/
3.kmerfreq gce

kmerfreq -p OUTPUT -t 10 reads_files.lib
-p 输出文件夹
-t 线程数
-k kmer值
将测序数据的路径写入reads_files.lib文件中

gce

less AF.kmer.freq.stat | grep "#Kmer indivdual number" 
less AF.kmer.freq.stat | perl -ne 'next if(/^#/ || /^\s/); print; ' | awk '{print $1"\t"$2}' > AF.kmer.freq.stat.2colum 

对于纯合和接近纯合的基因组

gce -g 173854609857 -f AF.freq.stat.2colum >gce.table 2>gce.log

对于杂合基因组

gce -g 173854609857 -f AF.freq.stat.2colum -c 75 -H 1 >gce2.table 2>gce2.log

查询植物基因组大小的网站：http://data.kew.org/cvalues
查询动物基因组大小的网站：http://www.genomesize.com

引用：
https://github.com/gmarcais/Jellyfish
https://kat.readthedocs.io/en/latest/using.html
https://github.com/fanagislab/GCE