组装结果统计

下面对组装得到fasta 格式基因组序列进行长度、N50 等统计。

使用assembly-stats对组装结果进行统计

下载地址：https://github.com/sanger-pathogens/assembly-stats

assembly-stats K41.scafSeq \ #输入组装结果文件
 > out.N50.stat #输出文件

cat out.N50.stat
stats for K41.scafSeq
sum = 4551705, n = 166, ave = 27419.91, largest = 236234 N50 = 95421, n = 15
N60 = 78576, n = 21
N70 = 57780, n = 28 N80 = 42471, n = 36 N90 = 29527, n = 49 N100 = 100, n = 166 N_count = 428
Gaps = 19

gc-depth 分析

得到组装结果后，可以通过将测序数据比对回拼接结果，将全基因组划分窗口，统计每个窗口的平均GC 含量和覆盖深度，得到gc-depth 图。

#构建基因组bwa index
bwa index genome.fasta

#比对并排序
bwa mem -t 4 \#线程数
genome.fasta ./ecoli_R1.fastq.gz ./ecoli_R2.fastq.gz | \ #输入文件
samtools sort - -o aln_sort.bam

#定义变量
genome=genome.fasta ## 基因组文件
bam=aln_sort.bam ## 比对结果文件
prefix=gcdep ## 输出结果前缀
window=500 ## 窗口大小
step=250 ## 步长

#计算基因组序列长度
seqtk comp $genome | awk '{print $1"\t"$2}' > $prefix.len

#划分窗口 生成bed文件
bedtools makewindows -w $window -s $step -g $prefix.len > $prefix.window.bed

#按窗口提取序列并计算gc含量
seqtk subseq $genome $prefix.window.bed > $prefix.window.fasta
seqtk comp $prefix.window.fasta |awk '{print $1 "\t" ($4+$5)/($3+$4+$5+$6) }' > $prefix.window.gc

#按窗口计算平均深度
bedtools coverage -b aln_sort.bam -a gcdep.window.bed -mean | awk '{print $1":"$2+1"-"$3"\t"$4}' > $prefix.window.depth

#绘图
Rscript run_gcdep.R $prefix.window.gc $prefix.window.depth $prefix.pdf 0 0.8 0 500