前言：

最近組裝了一種病原體的基因組，基因組大小為610kb客们，結(jié)果發(fā)現(xiàn)在300,000-400,000之間發(fā)現(xiàn)很多的Gap區(qū)域崇决，需要找一下原因材诽。因?yàn)槭怯枚鷶?shù)據(jù)測的，我先推測的原因是基因組這個(gè)區(qū)域有可能GC含量比較高恒傻，那下載一下它的基因組岳守，看一下，找到了bedtools工具碌冶，發(fā)現(xiàn)這個(gè)軟件功能十分強(qiáng)大，bedtools總共有二三十個(gè)工具/命令來處理基因組數(shù)據(jù)涝缝。比如：根據(jù)bed中的位置信息提取目標(biāo)基因及其上下游序列扑庞；統(tǒng)計(jì)基因組不同區(qū)間的GC含量；提取gff文件的所有基因位置,并轉(zhuǎn)換成bed格式拒逮，計(jì)算覆蓋度（coverage）（coverageBed罐氨，genomeCoverageBed）等等。

一．Bedtools評估基因組二代數(shù)據(jù)覆蓋度

分析思路：把這個(gè)基因組文件按照100,000bp大小滩援，按照堿基位置分割成6個(gè)文件栅隐，然后分別計(jì)算不同區(qū)間的GC含量。

1.軟件下載安裝：

https://github.com/arq5x/bedtools2

wget?https://github.com/arq5x/bedtools2/releases/download/v2.30.0/bedtools-2.30.0.tar.gz

tar?-zxf?bedtools-2.30.0.tar.gzcd?bedtools2/make

然后用help看一下玩徊，安裝成功租悄。

?2.軟件使用：

準(zhǔn)備基因組文件，如果是多個(gè)序列文件恩袱，比如我的基因組文件泣棋，先計(jì)算各個(gè)contig/scaffold的長度：

使用python里面的pyfaidx模塊的faidx命令，代碼如下：

conda activate py36 (我自己的一個(gè)py36小環(huán)境)

pip install pyfaidx

faidxminia_k81.contigs.fa -i chromsizes > size.genome

? ?結(jié)果如下：

劃分窗口：

/public/home/rp1016swf/rp1016swf/software/bedtools2/bin/bedtools makewindows -g sizes.genome -w 50000 > windows. Bed

-g sizes.genome是要?jiǎng)澐值幕蚪M畔塔，格式為兩列：染色體潭辈、染色體長度-w?50000為窗口大小為5wwindows.bed為輸出文件，格式為三列：染色體澈吨、區(qū)間開始位點(diǎn)把敢、區(qū)間結(jié)束位點(diǎn)。

?統(tǒng)計(jì)窗口內(nèi)的GC含量：

/bedtools2/bin/bedtools?nuc?-fi?ViralProj237323_genomic.fna?-bed?windows.bed?|?cut?-f?1-3,5?>?5w_gc.bed

統(tǒng)計(jì)窗口內(nèi)的平均覆蓋深度

bedtools?coverage?-a?windows.bed?-b?SRR081241.sorted.bam?>?RR081241.depth.txt

bedtools?coverage對劃分好的每個(gè)滑動(dòng)窗口進(jìn)行reads數(shù)（depth)的統(tǒng)計(jì)谅辣。

-a windows為上一步劃分好的區(qū)間

-SRR081241.sorted.bam為測序數(shù)據(jù)mapping到參考基因組的比對文件

HG00096.depth.txt為統(tǒng)計(jì)結(jié)果的輸出文件修赞，格式為7列：染色體、區(qū)間起始位點(diǎn)桑阶、區(qū)間結(jié)束位點(diǎn)榔组、該區(qū)間內(nèi)的reads數(shù)、該區(qū)間內(nèi)的堿基數(shù)联逻、區(qū)間大小搓扯、該區(qū)間的平均覆蓋度

生成的txt文件共有7列，分別為序列編號包归、起始位置锨推、結(jié)束位置、reads數(shù)、堿基數(shù)换可、區(qū)間大小椎椰、平均覆蓋深度

二、Bowtie2+Samtools評估基因組二代數(shù)據(jù)比對率

Bowtie下載安裝

wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.5.1/bowtie2-2.5.1-linux-x86_64.zip/

1.先是構(gòu)建索引：

#bowtie2-build -f ./minia_k81.contigs.fa -p ./minia_k81.contigs?-p索引文件前綴名

2.bowtie2比對及samtools轉(zhuǎn)為bam文件,并根據(jù)比對情況進(jìn)行提取

bwa比對生成的為sam（sequence Alignment mapping)文件沾鳄，將SAM轉(zhuǎn)換為二進(jìn)制對應(yīng)的BAM格式慨飘。二進(jìn)制格式對于計(jì)算機(jī)程序來說更容易使用。要將SAM轉(zhuǎn)換為BAM译荞，我們使用samtools view命令瓤的。

bowtie2 -x ./minia_k81.contigs -p 20 -1 20220652_mapped_P1.fq -2 20220652_mapped_P2.fq -S  ./minia_k81.contigs.sam

3.準(zhǔn)備序列比對后生成的 bam 文件或者 .sam 文件

samtools view -bS minia_k81.contigs.sam > minia_k81.contigs.bam

4.統(tǒng)計(jì)序列比對情況

samtools flagstat scaffolds.bam > flagstatstat.txt

本文使用 文章同步助手 同步