resequencing data analyses

4.png

1. mapping 率

#這里使用的是sort過(guò)后的bam；
#-@ 15 線程數(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

2. sequencing depth

#使用sort過(guò)后的bam；
#-a 輸出所有位點(diǎn)，包括零深度的位點(diǎn)；
#沒(méi)有找到線程數(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ū)別

bedtools 中g(shù)etfasta，根據(jù)坐標(biāo)區(qū)域來(lái)從基因組里面提取fasta序列
bedtools 用法大全

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