1,Fastq数据质控

2,Fastq转化成bam,包含头文件

bwa aln ref.fa test_1.fq > test_1.sai
bwa aln ref.fa test_2.fq > test_2.sai
bwa sampe ref.fa -r "@RG\tID:<ID>\tLB:<LIBRARY_NAME>\tSM:<SAMPLE_NAME>\tPL:ILLUMINA" test_1.sai test_2.sai test_1.fq test_2.fq > test.sam

3,sam 转化成bam,如果SAM文件中有header @SQ lines。

samtools view -b -S test.sam > test.bam
##如果没有header时: samtools faidx ref.fa
## samtools view -bt ref.fa.fai test.sam > test.bam

4,sort bam

samtools sort test.bam > test.sorted.bam

或者:

java -jar picard.jar SortSam I=test.bam O=test.sorted.bam SORT_ORDER=coordinate

5, 标记重复

java -jar picard.jar MarkDuplicate ....

6, index 一下

samtools index test.sorted.repeatmark.bam 

7,Base Quality Score Recalibration

....

8, 使用GATK检测SNP

java -jar GenomeAnalysisTK.jar glm SNP -R ref.fa -T UnifiedGenotyper -I test.sorted.repeatmark.bam -o test.raw.vcf

使用samtools和bcftools检测SNP

samtools faidx ref.fa
samtools mpileup -d 1000 -DSugf test.sorted.repeatmark.bam > test.raw.vcf ##(samtools mpileup -vf 。。。) bcftools view -Nvcg -d 1000 test.raw.vcf > test.snp.vcf ##(我的软件运行这步会出错,用下面两行代码代替)
bcftools call -mv test.raw.vcf > test.raw_varient.vcf
bcftools filter -s LowQual -e '%QUAL<20 || DP>100' test.raw_varient.vcf > test.filt_varient.vcf

##也有直接用perl 脚本实现。在使用bcftools 得到variant calling变异后的结果后。需要对结果再次进行过滤,主要依据对比结果中的第8列消息,其中的DP4最为重要,对应的提供了四列:1, 比对结果和正链一致的reads数;2, 比对结果和负链一致的reads数;3, 比对结果在正链的variant 上的reads数;4, 比对结果在负链的variant上的reads数。当设定(value3+value4)大于某一阈值,才算是variant。

bcftools检测生成的vcf格式有10列。1,参考序列名。2,variant所在的left-most位置。3,variant的ID,(默认未设置,用“.”表示)。4,参考序列的allele。5,variant的allele(有多个alleles,则用“,”分隔)。6,variant/reference Quality。7,FILTers applied。8,varient的信息,使用分号隔开。9,Format of the genotype fields, seperated by colon (optional)。10,Sample genotypes and per-sample information(optional)。

bcftools 的第8列中显示了对variants的信息描述,其中Tag的描述如下:

Tad  Format  Description

AF1  double  Max-likelihood estimate of the site allele Frequency (AF)of the first ALT allele

DP   int    Raw read depth (without quality filtering)

DP4  int[4]    high-quality reference forward base, ref reverse, alternate for and alt rev bases

FQ  int     consensus quality. Positive: sample genotypes different; negative: otherwise

MQ	int	Root-Mean-Square mapping quality of covering reads
PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
QCHI2 int Phred-scaled PCHI2
RP int # permutations yielding a smaller PCHI2
CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
UGT string Most probable genotype configuration without the trio constraint
CGT string Most probable configuration with the trio constraint

使用bcftools过滤掉不可靠的位点:

bcftools filter的参数:

-e -exclude 主要用于表达式方式去除匹配上的位点,这个参数很关键,过滤需要此表达式

-g -SnpGap 过滤INDEL附近的snp位点,比如-SnpGap 5 则过滤INDEL附近5个碱基距离内的SNP

-G -IndelGap 过滤INDEL附近的INDEL位点

-o -output 输出文件的名称

-O -output-type 输出的格式,一般z和v都行

-s -soft-filter 将过滤掉的位点用字符串注释

-S -set-GTs setgenotypes of failed samples to missing value (.) or reference allele (0) (将不符合要求的个体基因改为".")

eg:过滤QUAL小于10,DP值小于5,INDEL附近的位点

bcftools filter -O v  -o test.filter_variant.vcf -s LOWQUAL -e 'QUAL<10 || FMT/DP < 5' --SnpGap 5 --set-GTs .  test.vcf.gz

提取过滤后的SNP位点

bcftools view -v snps test.filter_variance.vcf > test.snp_filter.vcf

或者在vcf文件中的INFO列里,如果是INDEL的话,会标注出INDEL,因此提取SNP也可以:

grep -v 'INDEL'  test.filter_variance.vcf > test.snp_filter.vcf

注意:

|| 与 | 区别:都表示“或”运算, 但是 || 运算符第一个表达式成立的话,后面的表达式不运算,直接返回。而 |  对所有表达式都判断。 && 与 & 的区别同理

参考:https://www.cnblogs.com/xiaofeiIDO/p/6857745.html

   http://www.bioinfo-scrounger.com/archives/248

05-08 14:51