組裝結(jié)果統(tǒng)計(jì)
下面對(duì)組裝得到fasta 格式基因組序列進(jìn)行長(zhǎng)度、N50 等統(tǒng)計(jì)即纲。
使用assembly-stats對(duì)組裝結(jié)果進(jìn)行統(tǒng)計(jì)
下載地址：https://github.com/sanger-pathogens/assembly-stats
assembly-stats K41.scafSeq \ #輸入組裝結(jié)果文件
> 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 分析
得到組裝結(jié)果后啦撮，可以通過(guò)將測(cè)序數(shù)據(jù)比對(duì)回拼接結(jié)果谭网，將全基因組劃分窗口，統(tǒng)計(jì)每個(gè)窗口的平均GC 含量和覆蓋深度赃春，得到gc-depth 圖蜻底。
#構(gòu)建基因組bwa index
bwa index genome.fasta
#比對(duì)并排序
bwa mem -t 4 \#線(xiàn)程數(shù)
genome.fasta ./ecoli_R1.fastq.gz ./ecoli_R2.fastq.gz | \ #輸入文件
samtools sort - -o aln_sort.bam
#定義變量
genome=genome.fasta ## 基因組文件
bam=aln_sort.bam ## 比對(duì)結(jié)果文件
prefix=gcdep ## 輸出結(jié)果前綴
window=500 ## 窗口大小
step=250 ## 步長(zhǎng)
#計(jì)算基因組序列長(zhǎng)度
seqtk comp $genome | awk '{print $1"\t"$2}' > $prefix.len
#劃分窗口 生成bed文件
bedtools makewindows -w $window -s $step -g $prefix.len > $prefix.window.bed
#按窗口提取序列并計(jì)算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
#按窗口計(jì)算平均深度
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
