resequencing data analyses
WGS要分析什么
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
samtools markdup -@ 15 4AL_fm_sort.bam 4AL_fm_sort_markdup.bam
# 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 4AL_fm_sort_markdup.bam |bcftools call -vmO z -o 4AL_raw.vcf.gz
Pileup格式介绍
[samtools]mpileup命令简介
利用samtools mpileup和bcftools进行SNP calling
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
网友评论