resequencing data analyses
4.png
1. mapping 率
#这里使用的是sort过后的bam；
#-@ 15 线程数是15
nohup samtools flagstat -@ 15 4AL_reseq.sorted.bam >/home/huawei/raw_data/wu/WGS/1_sort_bam/4AL_reseq.sorted.bam.txt &
#查看生成内容
cat 4AL_reseq.sorted.bam.txt
1.png
samtools flagstat统计bam文件比对结果
2. sequencing depth
#使用sort过后的bam；
#-a 输出所有位点，包括零深度的位点；
#没有找到线程数；
nohup samtools depth 4AL_reseq.sorted.bam -a >/home/huawei/raw_data/wu/WGS/1_sort_bam/4AL_reseq.sorted.bam.depth &
第一列尾染色体号；第二列碱基；第三列为碱基上的测序深度
sequencing depth
参考：
samtools常用命令详解
3. read coverage
#-bga 在BEDGRAPH format中显示所有位置的genome coverage
#-ibam 输入文件是bam格式，必需经过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
#
解读结果：
`-bga` Reporting genome coverage for *all* positions in 在BEDGRAPH格式中显示”所有“位置的基因组覆盖度。-bg是实际覆盖在基因组上的区段，-a表示包括那些0 覆盖的区段。
read coverage结果
官方解读bedtools中genomecov的用法：
genomecov
官方配图：
coverage
2-3测序深度和覆盖度的区别：
![]()
测序深度和覆盖度的区别
测序深度和覆盖度的区别
bedtools 中getfasta，根据坐标区域来从基因组里面提取fasta序列
bedtools 用法大全
bedtools, vcftools, bcftools的功能区别
4. samtools mpileup +bcftools call 找SNP
#用mapping完后，samtools的sam转bam出的文件，默认按照name排序，标记,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进行排序
samtools sort -@ 15 -o 4AL_fm_sort.bam 4AL_fixmat.bam
#标记重复，但不要删除，samtools markdup [options] input output
#这个一直报错，报错如下，讲真我是run samtools fixmate，所以后面没有用samtools markdup，如picard或是gatk同样可以进行，
samtools markdup -@ 15 4AL_fm_sort.bam 4AL_fm_sort_markdup.bam
#用gatk来做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 来做markduplicate
java -jar picard -I 4AL_fm_sort.bam -M 4AL_fm_sort_mk_metrics.txt -O 4AL_fm_sort_mkdup.bam
报错
未完成7.24
# samtools mpileup 生成BCF、VCF文件或者pileup一个或多个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命令简介
利用samtools mpileup和bcftools进行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
