测序原始数据处理-质控

SAM文件详解：学习资源 http://davetang.org/wiki/tiki-index.php?page=SAM

      http://blog.csdn.net/u014182497/article/details/51691743

trim_galore --phred33 --gzip --length 20 --paired --trim-n -o ./ gDNA-Male_S2_L001_R1_001.fastq.gz gDNA-Male_S2_L001_R2_001.fastq.gz

fastqc gDNA-Male_S2_L001_R2_https://wenku.baidu.com/view/fb5c4dd4763231126fdb116b.html001_val_2.fq.gz -t 20 --noextract -o ./ *_val_1.fq.gz *_val_2.fq.gz

FastQc 相关学习

Fastqc 步骤后会得到html文件 

文件

***一般情况，Per base sequence quality和Per base GC content通过，测序质量基本合格，可以用于后续分析。

1.Per base sequence quality

quality就是Fred值，-10*log10(p)，p为测错的概率。所以一条reads某位置出错概率为0.01时，其quality就是20。图像如下面例子：

横轴代表位置，纵轴quality。红色表示中位数，黄色是25%-75%区间，触须是10%-90%区间，蓝线是平均数。

若任一位置的下四分位数低于10或中位数低于25，报"WARN"；若任一位置的下四分位数低于5或中位数低于20，报"FAIL".

2.Per Base Sequence Content

对所有reads的每一个位置，统计ATCG四种碱基（正常情况）的分布：

 

横轴为位置，纵轴为百分比。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。当所有位置的碱基比例一致的表现出bias时，即四条线平行但分开，往往代表文库有bias (建库过程或本身特点)，或者是测序中的系统误差。

当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"。

点击打开链接

点击打开链接

bwa mem -t 20 -R '@RG\tPL:Illumina\tID:foo\tSM:T' /home/ref/hg19/gatk/ucsc.hg19.fasta *_val_1.fq.gz *_val_2.fq.gz >pe.sam
###bwa 三种算法 bwa backtrack,bwa sw,baw men(最新的)。对于上述三种算法，首先需要使用索引命令构建参考基因组的索引，用于后面的比对。
#所以，使用BWA整个比对过程主要分为两步，第一步建索引，第二步使用BWA MEM进行比对。
# 建立索引：bwa index [-p prefix] [-a algoType] <ref.fa> 
# 参数详解：ref.fa——参加基因组文件，作为输入文件；
# p——输出文件前缀；
# a——构建索引的算法；包括两个算法，分别是is和bwtsw。对于参考基因组文件大于2G的使用bwtsw算法，使用bwtsw算法必须保证参考基因组文件大小大于10M。
#
# 比对：bwa mem -R '@RG\tID:foo\tSM:bar\tLB:library1' <ref.fa> <read1.fa> <read2.fa>  >  lane.sam
# 参数详解：
# R——设置reads标头，“\t”分割。 例如：’@RG\tID:foo\tSM:bar’；
# t——设置线程数；
# M——将较短的split hits标记为secondary，与picard兼容。
# f——文件输出  (>pe.sam)也可以的








java1.8 -jar $picard SortSam \
      I=pe.sam \
      O=sorted.bam \
      SORT_ORDER=coordinate


###对bam文件进行sort排序处理
#这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。可以使用picard-tools中SortSam完成。


#java -Xmx4g  设置jvm内存大小


#这一步对bam文件进行排序的时候用的是picard里面的SortSam这个工具，而不是直接用是samtools来进行sort，
#因为samtools sort排完序以后不会更新BAM文件的头文件，SO这个tag显示的可能还是unsorted，下游程序处理的时候可能会报错；SORT_ORDER这边设置成coordinate表示按照染色体位置进行排序






samtools stats sorted.bam >sorted.bam.stats
 -p bamstats/ sorted.bam.stats




java1.8 -jar $picard MarkDuplicates \
I=sorted.bam \
O=marked_duplicates.bam \
M=marked_dup_metrics.txt


#这一步是对每个sample的一些重复序列进行一些处理，这些重复的序列可能是因为PCR扩增的时候引入的一些引物序列，
#容易干扰下游结果，这个操作可以通过picard里面的MarkDuplicates工具实现，MarkDuplicates处理完成不会删除这些reads而是加上一个标签，
#带有这些标签的reads在后续处理的时候会被跳过。(注：dedup这一步只要在library层面上进行就可以了，例如一个sample如果建了多个库的话，
#对每个库进行dedup即可，不需要把所有库合成一个sample再进行dedup操作)


#用picard里面的MarkDuplicates工具去除PCR duplicates，这一步会生成
# 两个文件，一个dedup好的bam文件以及一个纯文本文件，这边起名叫
# marked_dup_metrics，里面包含duplicates的一些统计信息


java1.8 -jar $picard BuildBamIndex \
I=marked_duplicates.bam






java1.8 -jar /home/biosoftware/install_pkg/GenomeAnalysisTK.jar \
-T RealignerTargetCreator \               #用RealignerTargetCreator找到需要重新比对的区域,输出文件intervals
-R /home/ref/hg19/gatk/ucsc.hg19.fasta \  #参考基因组，参考基因组特别重要
-I marked_duplicates.bam \
-known /home/database/gatk/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-known /home/database/gatk/1000G_phase1.indels.hg19.sites.vcf \
-o dedup.realn.intervals 




java1.8 -jar /home/biosoftware/install_pkg/GenomeAnalysisTK.jar \
   -T IndelRealigner \                       ##用RealignerTargetCreator找到需要重新比对的区域,输出文件intervals
   -R /home/ref/hg19/gatk/ucsc.hg19.fasta \
   -I marked_duplicates.bam \
   -known /home/database/gatk/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
   -known /home/database/gatk/1000G_phase1.indels.hg19.sites.vcf \
   -targetIntervals dedup.realn.intervals \
   -o realigned.bam


java1.8 -jar /home/biosoftware/install_pkg/GenomeAnalysisTK.jar \
        -T BaseRecalibrator \
        -R /home/ref/hg19/gatk/ucsc.hg19.fasta \
        -I realigned.bam \
        -knownSites /home/database/gatk/dbsnp_138.hg19.vcf \
        -knownSites /home/database/gatk/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
        -o recal_data.grp


java1.8 -jar /home/biosoftware/install_pkg/GenomeAnalysisTK.jar \
-T PrintReads \
-R /home/ref/hg19/gatk/ucsc.hg19.fasta \
-I realigned.bam \
-BQSR recal_data.table \
-o recal_reads.bam






java1.7 -jar /home/biosoftware/install_pkg/mutect-src/mutect/target/mutect-1.1.7.jar \
-T MuTect \
-R /home/ref/hg19/gatk/ucsc.hg19.fasta \
--dbsnp /home/database/gatk/dbsnp_138.hg19.vcf \
--cosmic /home/wyy/drop/hg19_cosmic_v81.vcf \
  --input_file:normal recal_reads.bam \
--input_file:tumor realigned.bam \
--out call_stats.out \
--coverage_file coverage.wig.txt