|
阅读：5554回复：0
全外显子分析（WES）的一般步骤-3 QC篇
这步主要是判断一下NGS数据格式，以及做一些QC
############################################################# # Step 1. check the quality score encoding ############################################################# # S - Sanger Phred+33, raw reads typically (0, 40) # X - Solexa Solexa+64, raw reads typically (-5, 40) # I - Illumina 1.3+ Phred+64, raw reads typically (0, 40) # J - Illumina 1.5+ Phred+64, raw reads typically (3, 40) # with 0=unused, 1=unused, 2=Read Segment Quality Control Indicator # L - Illumina 1.8+ Phred+33, raw reads typically (0, 41) less $read_1 | head -n400 | awk '{if(NR%4==0) printf("%s",$0);}' | od -A n -t u1 -v | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}' less $read_2 | head -n400 | awk '{if(NR%4==0) printf("%s",$0);}' | od -A n -t u1 -v | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}' ############################################################# # Step 2. Quality Control ############################################################# mkdir -p $OUTPUT/QC #$fastqc -o $OUTPUT/QC -noextract -f fastq -q $read_1 $read_2 ############################################################# # Step 3. Detect and Remove Contaminations ############################################################# # if Overrepresented Sequences failed, skewer can be called #$skewer -x ADAPTER1 -y ADAPTER2 -q 3 -r 0 -d 0 -t $CPUs -l 25 -o $OUTPUT/OUTPUT_BASENAME $read_1 $read_2 ############################################################# # Step 4. Remove Low Quality Reads ############################################################# # run this only if it's needed #$sickle pe -f $read_1 -r $read_2 -t sanger -o $OUTPUT/read1_trim.fastq -p $OUTPUT/read2_trim.fastq -s $OUTPUT/single_trim.fastq -q 30 -l 25 |
|
