readme

README.md

@@ -29,7 +29,7 @@ choppy install renluyao/Quality_control
 
 ###1. 原始数据质量控制
 
-#### [Fastqc](<https://www.bioinformatics.babraham.ac.uk/projects/fastqc/>)
+#### [Fastqc](<https://www.bioinformatics.babraham.ac.uk/projects/fastqc/>) v0.11.5
 
 FastQC是一个常用的测序原始数据的质控软件，主要包括12个模块，具体请参考[Fastqc模块详情](<https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/>)。
 
@@ -37,7 +37,7 @@ FastQC是一个常用的测序原始数据的质控软件，主要包括12个模
 fastqc -t <threads> -o <output_directory> <fastq_file>
 ```
 
-#### [Fastq Screen](<https://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/>)
+#### [Fastq Screen](<https://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/>) 0.12.0
 
 Fastq Screen是检测测序原始数据中是否引⼊入其他物种，或是接头引物等污染，⽐比如，如果测序样本
 是⼈人类，我们期望99%以上的reads匹配到⼈人类基因组，10%左右的reads匹配到与⼈人类基因组同源性
@@ -54,7 +54,7 @@ fastq_screen --aligner <aligner> --conf <config_file> --top <number_of_reads> --
 
 ###2. 比对后数据质量控制
 
-#### [Qualimap](<http://qualimap.bioinfo.cipf.es/>)
+#### [Qualimap](<http://qualimap.bioinfo.cipf.es/>) 2.0.0
 
 Qualimap是一个比对指控软件，包含Picard的MarkDuplicates的结果和sentieon中metrics的质控结果。
 
@@ -64,7 +64,7 @@ qualimap bamqc -bam <bam_file> -outformat PDF:HTML -nt <threads> -outdir <output
 
 ###3. 突变检出数据质量控制
 
-#### [Hap.py](<https://github.com/Illumina/hap.py>)
+#### [Hap.py](<https://github.com/Illumina/hap.py>) v0.3.9
 
 hap.py是将被检测vcf结果与benchmarking对比，计算precision和recall的软件，它考虑了vcf中[突变表示形式的多样性](<https://genome.sph.umich.edu/wiki/Variant_Normalization>)，进行了归一化。
 
@@ -72,7 +72,7 @@ hap.py是将被检测vcf结果与benchmarking对比，计算precision和recall
 hap.py <truth_vcf> <query_vcf> -f <bed_file> --threads <threads> -o <output_filename>
 ```
 
-#### [Jaccard index](<https://en.wikipedia.org/wiki/Jaccard_index>)
+#### [Jaccard index](<https://en.wikipedia.org/wiki/Jaccard_index>) (rtg-tools 3.10.1)
 
 Jaccard index是两个集合的交集除以两个集合的并集。这个计算用的是rtg-tools的vcfeval，它将两个vcf中的突变表示格式进行了统一，即两个看似不同的SNV或者INDEL实际上是指同一个突变，这种情况经常发生在重复区域。
 
@@ -82,7 +82,7 @@ rtg vcfeval -b <one_vcf> -c <another_vcf> -o <output_directory> -t <sdf_file>
 
 `-t` sdf文件是rtg-tools专用的注释文件，由fasta文件转换来
 
-#### [VCF statistics](<https://github.com/RealTimeGenomics/rtg-tools>)
+#### [VCF statistics](<https://github.com/RealTimeGenomics/rtg-tools>) (rtg-tools 3.10.1)
 
 该计算用的是rtg-tools的vcfstats，对vcf中的SNV和INDEL的数目进行了统计。
 
@@ -90,7 +90,18 @@ rtg vcfeval -b <one_vcf> -c <another_vcf> -o <output_directory> -t <sdf_file>
 rtg vcfstats <vcf_file>
 ```
 
+### 4. 质控数据整合
+
+####[MultiQC](<http://multiqc.info/>) v1.8
+
+以上的质控软件的输出都是单个文件的质控结果，multiqc可以将这些结果进行汇总，以及网页可视化展示
+
+```bash
+multiqc <input_directory>
+```
+
 ## App输入变量与输入文件
+
 在安装了APP之后，输入一下命令得到需要准备的文件
 
 ```bash
@@ -135,53 +146,93 @@ oss://chinese-quartet/quartet-test-data/fastqfiles/Fudan_DNA_LCL5_R1.fastq.gz	os
 **sample_id**是choppy app内置的识别标志，写1即可
 
 ## App输出文件
+以下task的输出都只包含了单个文件的质控结果
+
+**fastq.wdl**
+
+- _fastqc.html
+- _fastqc.zip
+
+**fastqscreen.wdl**
+
+- _screen.png
+- _screen.txt
+- _screen.html
+
+**bamqc.wdl**
+
+- _qualimap.zip
+
+**benchmark.wdl**
 
+- .rtg.vcf.gz
+- .rtg.vcf.gz.tbi
+- .vcf.gz
+- .vcf.gz.tbi
+- .roc.all.csv.gz
+- .roc.Locations.INDEL.csv.gz
+- .roc.Locations.INDEL.PASS.csv.gz
+- .roc.Locations.SNP.csv.gz
+- .roc.Locations.SNP.PASS.csv.gz
+- .summary.csv
+- .extended.csv
+- .metrics.json.gz
+
+**jaccard_index.wdl**
+
+- summary.txt
+
+**vcfstat.wdl**
+
+- onestats.txt
+
+主要查看以下三个task的结果，汇总了所有样本的结果
+
+**mergeJI.wdl**
+
+- result.txt
+
+**mergeNum.wdl**
+
+- vcfstats.txt
+
+**multiqc**
+
+- multiqc_report.html
+- multiqc.log
+- multiqc_data.json
+- multiqc_fastq_screen.txt
+- multiqc_fastqc.txt
+- multiqc_general_stats.txt
+- multiqc_happy_data.json
+- multiqc_sources.txt
 
 ## 结果展示与解读
-GSEA结果解读示例：
-
-> ### 1. Enrichment score（ES）
->
-> ES是GSEA最初的结果，反应全部杂交data排序后，在此序列top或bottom富集的程度。
-> ES原理：扫描排序序列，当出现一个功能集中的gene时，增加ES值，反之减少ES值，所以ES是个动态值。最终ES的确定是讲杂交数据排序序列所在位置定义为0，ES值定义为距离排序序列的最大偏差.
-> - ES为正，表示某一功能gene集富集在排序序列前方
-> - ES为负，表示某一功能gene集富集在排序序列后方。
-> 图中的最高点为此通路的ES值，中间表示杂交数据的排序序列。竖线表示此通路中出现的芯片数据集中的gene。
->
-> ### 2. NES
-> 
-> 由于ES是根据分析的数据集中的gene是否在一个功能gene set中出现来计算的，但各个功能gene set中包含的gene数目不同，且不同功能gene set与data之间的相关性也不同，因此，比较data set在不同功能gene set中的富集程度要对ES进行标准化处理，也就是NES
-> NES=某一功能gene set的ES/数据集所有随机组合得到的ES平均值
-> NES是主要的统计量。
-> 
-> ### 3. FDR
-> 
-> NES确定后，判断其中可能包含的错误阳性发现率。FDR=25%意味着对此NES的确定，4次可能错  1次。GSEA结果中，高亮显示FDR<25%的富集set。因为从这些功能gene中最可能产生有意义的假设，促进进一步研究。大多数情况下，选FDR<25%是合适的，但是，假如分析的芯片data set较少，选择的是探针随机组合而不是表型组合，若p不严格，那么应该选FDR<5%。一般而言，NES绝对值越大，FDR值就越小，说明富集程度高，结果可靠。
-> 
-> ### 4. 名义p值 nominal p-value
-> 
-> 描述的是针对某一功能gene子集得到的富集得分的统计显著性，显然，p越小，富集性越好。
-> 
-> **以上4个参数中，只有FDR进行了功能gene子集大小和多重假设检验矫正，而p值没有，因此，如果结果中有一个高度富集的功能gene子集，而其有很小的名义p-value和大的FDR意味着富集并不显著。**
-> 
-> 我的一个具体结果解读：
-> 
-> > 92/681 gene sets are  upregulated in PH
-> > 0 gene sets are significantly enriched at FDR<25%
-> > 1 gene sets are significantly enriched at n p-value <1%
-> > 1 gene sets are significantly enriched at n p-value <5%
-> 
-> 在选择的BP中，有681个gene sets，92个PH中上调，其中75%的正确率支持0条子集上调，1个BP的gene表达上调名义p值<0.01。总体结果并不理想。
-> 
-> ### 5. 备注
-> 
-> #### GSEA富集结果太少说明：
-> 
-> 无gene set被富集。可能是因为分析的样本太少，关注的生物信息太微弱，或正在分析的功能集不能很好代表你所关心的生物过程，但仍然可以看下top ranked gene sets，这些信息可能会为你的假说提供微弱的证据。当然也可以尝试考虑分析其他gene sets，或增加samples
-> 
-> #### GSEA富集结果太多说明：
-> 
-> 太多的功能子集被富集了。可能是因为很多的gene sets代表同一生物信号，这可以在gene sets中查看leading edge sbusets来查看。或者也可以查看具体区别进行加工，比如samples来自不同labs，操作者不一样等。
+
+#### 1. result.txt
+
+```bash
+#vcf1	#vcf2	#vcf1_vcf2	#True-pos-call-number	#False-pos-number	#False-neg-number
+oss://choppy-cromwell-result/test-choppy/wgs_quartettest_renluyao_0827/72f269f2-91b7-4fbe-bde7-99b2e1e3091c/call-Haplotyper/Fudan_DNA_LCL7_hc.vcf	oss://choppy-cromwell-result/test-choppy/wgs_quartettest_renluyao_0827/7a72d0e6-302d-43ca-b6b0-daeaa0236d06/call-Haplotyper/Fudan_DNA_LCL5_hc.vcf	LCL7_LCL5	4891550	11116	11116
+```
+
+在获得result.txt之后，jaccard index的计算公式如下
+$$
+JaccardIndex=\frac{TP}{TP+FP+FN}
+$$
+
+#### 2. vcfstats.txt
+
+```bash
+File	Failed Filters	Passed Filters	SNPs	MNPs	Insertions	Deletions	Indels	Same as reference	SNP Transitions/Transversions	Total Het/Hom ratio	SNP Het/Hom ratio	MNP Het/Hom ratio	Insertion Het/Hom ratio	Deletion Het/Hom ratio	Indel Het/Hom ratio	Insertion/Deletion ratio	Indel/SNP+MNP ratio
+ /cromwell_inputs/choppy-cromwell-result/test-choppy/qc_test_renluyao_0831/79830609-9bd2-4e0c-9483-0c3369052a9d/call-benchmark/shard-0/Fudan_DNA_LCL5_hc.rtg.vcf.gz	 0	 4904087	 4024591	 0	 421154	 445839	 12503	 0	 1.98 (3771780/1907484)	 1.44 (2897278/2006809)	 1.43 (2371579/1653012)	 - (0/0)	 1.37 (243816/177338)	 1.53 (269380/176459)	 - (12503/0)	 0.94 (421154/445839)	 0.22 (879496/4024591)
+```
+
+#### 3. MultiQC输出的结果
+
+下载之后将文件名glob**改成multiqc_data，即可打开multiqc_report.html查看可视化的结果，在multiqc_data中的整合结果txt文件，可以用于报告系统的输入。
+
+对应的fastqc、fastqscreen、qualimap、hap.py的结果解释请查询对应的官网。
 
 ## CHANGELOG
 **Version 1.0 - Auguest 30, 2019**
@@ -199,5 +250,13 @@ RNAseq和甲基化的质控流程待完善
 
 在Version 1.0中暂时还没有考虑没有技术重复的问题，可输入姐妹、父母、父女、母女的配对，计算同卵双胞胎、亲属关系和陌生人之间基因突变位点的一致性。
 
+**3. 怎么对该APP的输出结果进行可视化？**
+
+正在努力开发中
+
+**4. bam文件和vcf文件怎么获得？**
+
+在进行质控分析前，请先用标准化流程进行测序数据分析，详情查看choppy APP [huangyechao/wgs-germline](<http://choppy.3steps.cn/huangyechao/wgs-germline>)
+
 
 

tasks/bamqc.wdl

@@ -23,11 +23,5 @@ task bamqc {
 
 	output {
 		File zip = "${bamname}_qualimap.zip"
-		File genome_result = "result/genome_results.txt"
-		File pdf = "result/report.pdf"
-		File coverage = "result/raw_data_qualimapReport/coverage_histogram.txt"
-		File insert_size = "result/raw_data_qualimapReport/insert_size_histogram.txt"
-		File genome_fraction = "result/raw_data_qualimapReport/genome_fraction_coverage.txt"
-		File gc_dist = "result/raw_data_qualimapReport/mapped_reads_gc-content_distribution.txt"
 	}
 }