resequencing data analyses
4.png
1. mapping 率
#這里使用的是sort過(guò)后的bam矩父；
#-@ 15 線(xiàn)程數(shù)是15
nohup samtools flagstat -@ 15 4AL_reseq.sorted.bam >/home/huawei/raw_data/wu/WGS/1_sort_bam/4AL_reseq.sorted.bam.txt &
#查看生成內(nèi)容
cat 4AL_reseq.sorted.bam.txt
1.png
samtools flagstat統(tǒng)計(jì)bam文件比對(duì)結(jié)果
2. sequencing depth
#使用sort過(guò)后的bam锉桑；
#-a 輸出所有位點(diǎn)，包括零深度的位點(diǎn)窍株；
#沒(méi)有找到線(xiàn)程數(shù)民轴；
nohup samtools depth 4AL_reseq.sorted.bam -a >/home/huawei/raw_data/wu/WGS/1_sort_bam/4AL_reseq.sorted.bam.depth &
第一列尾染色體號(hào)；第二列堿基球订；第三列為堿基上的測(cè)序深度
sequencing depth
參考：
samtools常用命令詳解
3. read coverage
#-bga 在BEDGRAPH format中顯示所有位置的genome coverage
#-ibam 輸入文件是bam格式后裸，必需經(jīng)過(guò)sort
nohup genomeCoverageBed -ibam 4AL_reseq.sorted.bam -g /home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0 -bga >/home/huawei/raw_data/wu/WGS/1_sort_bam/4AL_cov.bedgraph &
cat 4AL_cov.bedgraph|head -10
#
解讀結(jié)果：
`-bga` Reporting genome coverage for *all* positions in 在BEDGRAPH格式中顯示”所有“位置的基因組覆蓋度。-bg是實(shí)際覆蓋在基因組上的區(qū)段冒滩，-a表示包括那些0 覆蓋的區(qū)段微驶。
read coverage結(jié)果
官方解讀bedtools中g(shù)enomecov的用法：
genomecov
官方配圖：
coverage
2-3測(cè)序深度和覆蓋度的區(qū)別：
![]()
測(cè)序深度和覆蓋度的區(qū)別
測(cè)序深度和覆蓋度的區(qū)別
bedtools 中g(shù)etfasta，根據(jù)坐標(biāo)區(qū)域來(lái)從基因組里面提取fasta序列
bedtools 用法大全
bedtools, vcftools, bcftools的功能區(qū)別
4. samtools mpileup +bcftools call 找SNP
#用mapping完后，samtools的sam轉(zhuǎn)bam出的文件媒咳，默認(rèn)按照name排序瓦糕，標(biāo)記,samtools fixmate [options] input output
samtools fixmate -@ 15 /home/huawei/raw_data/YSQ/4AL_resequence/output/bam/4AL_reseq.bam 4AL_fixmat.bam
#將name排序的bam按照position進(jìn)行排序
samtools sort -@ 15 -o 4AL_fm_sort.bam 4AL_fixmat.bam
#標(biāo)記重復(fù)，但不要?jiǎng)h除款筑，samtools markdup [options] input output
#這個(gè)一直報(bào)錯(cuò)智蝠，報(bào)錯(cuò)如下，講真我是run samtools fixmate奈梳，所以后面沒(méi)有用samtools markdup杈湾，如picard或是gatk同樣可以進(jìn)行，
samtools markdup -@ 15 4AL_fm_sort.bam 4AL_fm_sort_markdup.bam
#用gatk來(lái)做markduplicate
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" MarkDuplicates -I 4AL_fm_sort.bam -M 4AL_fm_sort_mk_metrics.txt -O 4AL_fm_sort_mkdup.bam
#用picard 來(lái)做markduplicate
java -jar picard -I 4AL_fm_sort.bam -M 4AL_fm_sort_mk_metrics.txt -O 4AL_fm_sort_mkdup.bam
報(bào)錯(cuò)
未完成7.24
# samtools mpileup 生成BCF颈嚼、VCF文件或者pileup一個(gè)或多個(gè)bam文件
# bcftools call
samtools mpileup -ugf /home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta -t DP -t SP ../bam/4AL_fm_sort_mkdup.bam |bcftools call -vmO z -o 4AL_raw.vcf.gz
Pileup格式介紹
[samtools]mpileup命令簡(jiǎn)介
利用samtools mpileup和bcftools進(jìn)行SNP calling
Question: Why GATK and bcftools SNP calling different?
5. VCFTOOLs
vcftools --vcf X.vcf --remove-indels --out X.snps --recode --recode-INFO-all
vcftools --vcf X.vcf --keep-only-indels --out X.indel --recode --recode-INFO-all
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
bam="/home/huawei/raw_data/YSQ/4AL_resequence/output/bam/4AL_reseq.sorted.unique.markdup.add.bam
gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" HaplotypeCaller -R $ref -I $bam -O /home/huawei/raw_data/wu/WGS/VCF/4AL_resequence.vcf 1>${id}_log.HC 2>&1
