1.jellyfish

安裝jellyfish
$ wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.10.tar.gz
$ tar zxvf jellyfish-1.1.10.tar.gz
$ mkdir jellyfish
$ cd jellyfish-1.1.10
$ ./configure --prefix=$PWD/../jellyfish
如果安裝在當前目錄中，會報錯。
$ make -j 8
$ make install

運行jellyfish

cp denovo/newBGIseq500* .
gunzip  newBGIseq500_1.fq.gz
gunzip  newBGIseq500_2.fq.gz
#第一種方法
program/jellyfish count -m 17 -o kmer17  -s  268435456 -t 3 -c 8 -C newBGIseq500_1.fq newBGIseq500_2.fq
program/jellyfish stats -o kmer17.stats  kmer17_0
program/jellyfish histo -t 3  kmer17_0  > kmer17.histo
less kmer17.histo | awk '{sum +=$2*$1}END{print sum/41}'
#6.55551e+06
cat kmer17.stats kmer17.histo |perl -lane '{if(/Distinct:\s+(\d+)/){$total=$1/100}elsif(/\d+\s+(\d+)/){print "$1\t".($1/$total);}}' >kmer_freq_distribution
#畫圖
Rscript kmer_plot.R
`
png("kmer_distibution.png")
a <- read.delim("kmer_freq_distribution",head=F)
l<-seq(1,nrow(a),by=1 )
plot(x=l,y=a[,2],type ="l",col ="red",lwd=2,xlim=c(0,200),ylim=c(0,1.5),main="kmer頻率分布圖",xlab="depth",ylab="Frequency(%)")
#text(40,0.05,"38")
dev.off()
`



#第二種方法
program/jellyfish count -m 17 -o 17k1  -s  268435456 -t 3 -c 8 -C newBGIseq500_1.fq
program/jellyfish count -m 17 -o 17k2  -s  268435456 -t 3 -c 8 -C newBGIseq500_2.fq
program/jellyfish merge -o 17kk.jf 17k*

diff 17kk.jf kmer17_0
#2個文件一樣，沒有區(qū)別

image.png

2.使用 GCE 進行基因組大小評估

gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash \
 -k 17 -l read_list -i 450000000 -p species &> kmer_freq.log
gce-1.0.0/gce -f species.freq.stat -c 41 -H 1 -g 268799832 -m 1 -D 8 > species.table 2> species.log

下載GCE
wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
tar zxf gce-1.0.0.tar.gz -C /opt/biosoft

GCE 軟件包中主要包含 kmer_freq_hash 和 gce 兩支程序。前者用于進行 kmer 的頻數(shù)統(tǒng)計，后者在前者的結(jié)果上進行基因組大小的準確估算。

kmer_freq_hash 的常用參數(shù)：

Usage:  kmer_freq_hash [options]
    -k <int>    k-mer size(9~27), default=17
    -l <string> set read file list
    -c <double> set min accuracy rate of k-mer, set between 0~0.99, or -1 for unrestrained, default=-1
    -q <int>    set Quality value ascii shift, generally 64 or 33, default 64
    -r <int>    read length used to get k-mers, default=read's real length
    -a <int>    ignored length of the beginning of a read, default=0
    -d <int>    ignored length of the end of a read, default=0
    -g <int>    total bases used to get k-mers, default=all input bases
    -i <int>    initial size of hash table, default=1048576
    -p <string> set output prefix, default=output
    -o <int>    set whether ouput k-mer frequence file, 0 for no, 1 for yes, default=1
    -t <int>    thread number, default=1
    -L <int>    maximum read length, default=100
    -h      output help information to screen

Example: 
kmer_freq_hash -k 17 -i 450000000 -l fq.lst  2>kmerfreq.log;
kmer_freq_hash -k 19 -L 150 -i 600000000 -c 0.9 -t 8 -l fq.list  2>kmerfreq.log;

Attension: Please don't set -d and -r at the same time.

-k <17>
    設(shè)置 kmer 的大小。該值為 9~27，默認值為 17 。
-l string
    list文本文件，其中每行為一個fastq文件的路徑。
-t int
    使用的線程數(shù)，默認為 1 。
-i int
    初始的 hash 表大小，默認為 1048576。該值最好設(shè)置為 （kmer 的種類數(shù) / 0.75）/ 線程數(shù)。如果基因組大小為 100M，測序了 40M 個 reads，reads 的長度為 100bp，測序錯誤率為 1%，kmer的大小為 21，則kmer的種類數(shù)為100M+40M*100*1%*21=940M，若使用24線程，則該參數(shù)設(shè)置為 i=940M/0.75/24=52222222。
-p string
    設(shè)置輸出文件的前綴。
-o int
    是否輸出 k-mer 序列。1: yes, 0: no，默認為 1 。推薦選 0 以節(jié)約運行時間。
-q int
    設(shè)置fastq文件的phred格式，默認為 64。該值可以為 33 或 63。
-c double
    設(shè)置k-mer最小的精度，該值位于 0~0.99，或為 -1。 -1 表示不對 kmer進行過濾。設(shè)置較高的精度，可以用于過濾低質(zhì)量 kmer。精度是由 phred 格式的堿基質(zhì)量計算得來的。
-r int
    設(shè)置獲取 k-mer 使用到的 reads 長度。默認使用 reads 的全長。
-a int
    忽略read前面該長度的堿基。
-d int
    忽略read后面該長度的堿基。
-g int
    設(shè)置使用該數(shù)目的堿基來獲取 k-mers，默認是使用所有的堿基來獲取 k-mer。

運行kmer_freq_hash：

gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash \
 -k 17 -l read_list -i 450000000 -p species &> kmer_freq.log

kmer_freq_hash 的主要結(jié)果文件為 species.freq.stat。該文件有 2 列：第1列是kmer重復(fù)的次數(shù)，第二列是kmer的種類數(shù)。該文件有255行，第225行表示kmer重復(fù)次數(shù)>=255的kmer的總的種類數(shù)。該文件作為 gce 的輸入文件。
kmer_freq_hash 的輸出到屏幕上的信息結(jié)果保存到文件 kmer_freq.log 文件中。該文件中有粗略估計基因組的大小。其中的 Kmer_individual_num 數(shù)據(jù)作為 gce 的輸入?yún)?shù)。

less -S kmer_freq.log
**********Summary Table**********
Kmer_size       Kmer_individual_num     Kmer_species_num        Filter_kmer_num Kmer_depth(coverage)    Genome_size(rough)      Base_num        Base_depth(coverage)
17      268799832       25667758        0       41.2466 5993008 319999800       53.3955 3199998 100
#看不清
Kmer_size       17
Kmer_individual_num     268799832
Kmer_species_num        25667758
Filter_kmer_num             0
Kmer_depth(coverage)    41.2466
Genome_size(rough)      5993008
Base_num                     319999800
Base_depth(coverage)   53.3955 3199998 100

$ head species.freq.stat 
1   18626187
2   1256204
3   125673
4   28077
5   11014
6   6299
7   5330
8   5810
9   6887
10  7796                          

自己計算下Genome_size：peak=41
$ less -S species.freq.stat | awk '{sum += $1*$2} END {print sum/41}'
6.4344e+06   #這里是6.4M 比估算的6M要大

去掉第一行（1   18626187），計算為5980098,跟估算結(jié)果差不多
>a <-read.table("species.freq.stat")
>b <-a[,1]*a[,2]
> sum(b[2:length(b)])/41
[1] 5980098

gce 的使用：

 gce-1.0.0/gce -f species.freq.stat \
 -c 41 -g 268799832  -m 1 -D 8 > species.table 2> species.log

參數(shù)說明：

-f string
    kmer depth frequency file
-c int
    kmer depth frequency 的主峰對應(yīng)的 depth。gce 會在該值附近找主峰。
-g int
    總共的 kmer 數(shù)。一定要設(shè)定該值，否則 gce 會直接使用 -f 指定的文件計算 kmer 的總數(shù)。由于默認下該文件中最大的 depth 為 255，因此，軟件自己計算的值比真實的值偏小。同時注意該值包含低覆蓋度的 kmer。
-M int
    支持最大的 depth 值，默認為 256 。
-m int
    估算模型的選擇，離散型（0），連續(xù)型（1）。默認為 0，對真實數(shù)據(jù)推薦選擇 1 。
-D int
    precision of expect value，默認為 1。如果選擇了 -m 1，推薦設(shè)置該值為 8。
-H int
    使用雜合模式（1），不使用雜合模式（0）。默認值為 0。只有明顯存在雜合峰的時候，才選擇該值為 1 。
-b int
    數(shù)據(jù)是（1）否（0）有 bias。當 K > 19時，需要設(shè)置 -b 1 。

gce 的結(jié)果文件為 species.table 和 species.log 。species.log 文件中的主要內(nèi)容：

raw_peak        now_node        low_kmer        now_kmer        cvg     genome_size     a[1]    b[1]
41      5614870 21725832        242490544       40.8357 6.05044e+06     0.944654        0.844481

raw_peak： 覆蓋度為 41 的 kmer 的種類數(shù)最多，為主峰。
now_node： kmer的種類數(shù)。
low_kmer： 低覆蓋度的 kmer 數(shù)。
now_kmer： 去除低覆蓋度的 kmer 數(shù)，此值 = （-g 參數(shù)指定的總 kmer 數(shù)） - low_kmer 。
cvg：估算出的平均覆蓋度
genome_size：基因組大小，該值 = now_kmer / cvg 。
a[1]： 在基因組上僅出現(xiàn) 1 次的 kmer 之 種類數(shù)比例。
b[1]： 在基因組上僅出現(xiàn) 1 次的 kmer 之 數(shù)量比例。該值代表著基因組上拷貝數(shù)為 1 的序列比例。

如果使用 -H 1 參數(shù)，則會得額外得到如下信息：

for hybrid: a[1/2]=0.109694 a1=0.824146
kmer-species heterozygous ratio is about 0.0580297
for hybrid: b[1/2]=0.0680686 b1=0.77641
kmer-individual heterozygous ratio is about 0.0352335

Final estimation table:
raw_peak        now_node        low_kmer        now_kmer        cvg     genome_size     a[1/2]  a[1]    b[1/2]  b[1]
40      5621531 21685380        242491591       40.8306 6.05219e+06     0.109694        0.824146        0.0680686       0.77641

上面結(jié)果中，0.0580297 是由 a[1/2] 計算出來的。 0.0580297= a[1/2] / （ 2- a[1/2] ) 。
a[1/2]=0.109694 表示在所有的 uniqe kmer 中，有 0.109694 比例的 kmer 屬于雜合 kmer 。

此外，有 a[1/2] 和 b[1/2] 的值在最后的統(tǒng)計結(jié)果中。重復(fù)序列的含量 = 1 - b[1/2] - b[1] 。

則雜合率 = 0.0580297 / kmer_size 。 若計算出的雜合率低于 0.2%，個人認為測序數(shù)據(jù)應(yīng)該是純合的。這時候，應(yīng)該不使用 -H 1 參數(shù)。使用 -H 1 參數(shù)會對基因組的大小和重復(fù)序列含量估算造成影響。

參考：https://www.plob.org/article/9388.html

3.KmerFreq_AR計算基因組大小

 /home/zhaosl/biosoft/Assemblathon1_pipeline/KmerFreq_AR_v2.0 \
 -k 17 -t 20 -p test read_list

ls test*
test.freq.cz  test.freq.cz.len  test.freq.stat  test.genome_estimate

less -S test.genome_estimate
kmer    kmer_num        kmer_depth      genome_size     base_num        read_num        base_depth(X)   average_read_len        unique_kmer_number
17      268799832       41.2466 5993008 319999800       3199998 53.3955 100     25667758

kmer    17
kmer_num        268799832
kmer_depth      41.2466
genome_size     5993008
base_num        319999800
read_num        3199998
base_depth(X)   53.3955
average_read_len      100  
unique_kmer_number  25667758