split gvcf

README.md

+# WGS-germline Small Variants Quality Control Pipeline
+
+> Author： Run Luyao
+>
+> E-mail：18110700050@fudan.edu.cn
+>
+> Git: http://choppy.3steps.cn/renluyao/WGS_germline_datapotal.git
+>
+> Last Updates: 2020/0705
+
+## 安装指南
+
+```
+# 激活choppy环境
+source activate choppy
+# 安装app
+choppy install renluyao/WGS_germline_datapotal
+```
+
+## App概述——中华家系1号标准物质介绍
+
+建立高通量全基因组测序的生物计量和质量控制关键技术体系，是保障测序数据跨技术平台、跨实验室可比较、相关研究结果可重复、数据可共享的重要关键共性技术。建立国家基因组标准物质和基准数据集，突破基因组学的生物计量技术，是将测序技术转化成临床应用的重要环节与必经之路，目前国际上尚属空白。中国计量科学研究院与复旦大学、复旦大学泰州健康科学研究院共同研制了人源中华家系1号基因组标准物质（**Quartet，一套4个样本，编号分别为LCL5，LCL6，LCL7，LCL8，其中LCL5和LCL6为同卵双胞胎女儿，LCL7为父亲，LCL8为母亲**），以及相应的全基因组测序序列基准数据集（“量值”），为衡量基因序列检测准确与否提供一把“标尺”，成为保障基因测序数据可靠性的国家基准。人源中华家系1号基因组标准物质来源于泰州队列同卵双生双胞胎家庭，从遗传结构上体现了我国南北交界的人群结构特征，同时家系的设计也为“量值”的确定提供了遗传学依据。
+
+中华家系1号DNA标准物质的Small Variants标称值包括高置信单核苷酸变异信息、高置信短插入缺失变异信息和高置信参考基因组区。该系列标准物质可以用于评估基因组测序的性能，包括全基因组测序、全外显子测序、靶向测序，如基因捕获测序；还可用于评估测序过程和数据分析过程中对SNV和InDel检出的真阳性、假阳性、真阴性和假阴性水平，为基因组测序技术平台、实验室、相关产品的质量控制与性能验证提供标准物质和标准数据。此外，我们还可根据中华家系1号的生物遗传关系计算同卵双胞胎检测突变的一致性和符合四口之家遗传规律的一致率估计测序错误的比例，评估数据产生和分析的质量好坏。
+
+![](/Users/luyaoren/Documents/markdown_pictures/Quality_control/pictures/Picture1.png)
+
+该Quality_control APP用于全基因组测序（whole-genome sequencing，WGS）数据的质量评估，包括原始数据质控、比对数据质控和突变检出数据质控。
+
+## 流程与参数
+
+
+

tasks/benchmark.wdl

 	File vcf
 	File benchmarking_dir
 	File ref_dir
-	String sample = basename(vcf,".vcf")
+	String sample = basename(vcf,".splited.vcf")
 	String fasta
 	String docker
 	String cluster_config
 		/opt/rtg-tools/dist/rtg-tools-3.10.1-4d58ead/rtg index -f vcf ${sample}.rtg.vcf.gz
 
 		if [[ ${sample} =~ "LCL5" ]];then
-			/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL5.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL5.high.confidence.bed --threads $nt -o ${sample}
+			/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL5.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL5.high.confidence.bed --threads $nt -o ${sample} -r ${ref_dir}/${fasta}
 	    elif [[ ${sample} =~ "LCL6" ]]; then
-	    	/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL6.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL6.high.confidence.bed --threads $nt -o ${sample}
+	    	/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL6.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL6.high.confidence.bed --threads $nt -o ${sample} -r ${ref_dir}/${fasta}
         elif [[ ${sample} =~ "LCL7" ]]; then
-        	/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL7.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL7.high.confidence.bed --threads $nt -o ${sample}
+        	/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL7.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL7.high.confidence.bed --threads $nt -o ${sample} -r ${ref_dir}/${fasta}
 	    elif [[ ${sample} =~ "LCL8" ]]; then
-			/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL8.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL8.high.confidence.bed --threads $nt -o ${sample}
+			/opt/hap.py/bin/hap.py ${benchmarking_dir}/LCL8.afterfilterdiffbed.vcf.gz ${sample}.rtg.vcf.gz -f ${benchmarking_dir}/LCL8.high.confidence.bed --threads $nt -o ${sample} -r ${ref_dir}/${fasta}
         else
         	echo "only for quartet samples"
         fi		

tasks/fastqscreen.wdl

 	File read2
 	File screen_ref_dir
 	File fastq_screen_conf
-	String read1name = basename(read1,".fq.gz")
-	String read2name = basename(read2,".fq.gz")
+	String read1name = basename(read1,"\\.(fastq|fq)\\.gz$")
+	String read2name = basename(read2,"\\.(fastq|fq)\\.gz$")
 	String docker
 	String cluster_config
 	String disk_size

tasks/split_gvcf_files.wdl

 
 		sed -i '1,9d' name
 
-		for i in $(seq 1 $ncol);  do  cat ${gvcf}| cut -f1-9,$i > $i.splited.vcf; done
+		cat ${gvcf} | grep '#' > header
+		cat ${gvcf} | grep -v '#' > body
+		cat body | grep -w '^chr1\|^chr2\|^chr3\|^chr4\|^chr5\|^chr6\|^chr7\|^chr8\|^chr9\|^chr10\|^chr11\|^chr12\|^chr13\|^chr14\|^chr15\|^chr16\|^chr17\|^chr18\|^chr19\|^chr20\|^chr21\|^chr22\|^chrX' > body.filtered
+		cat header body.filtered > ${project}.filtered.g.vcf
+
+
+		for i in $(seq 10 $ncol);  do  cat ${project}.filtered.g.vcf | cut -f1-9,$i > $i.splited.vcf; done
 
 		ls *splited.vcf | sort -n | paste - name > rename
 
 		cat rename | while read a b
 		do
-			mv $a $b.vcf
+			mv $a $b.splited.vcf
 			if [[ $b.vcf =~ "LCL5_1" ]];then
 				cp $b.vcf ${project}.LCL5_1.vcf
 		    elif [[ $b.vcf =~ "LCL5_2" ]]; then
 		dataDisk: "cloud_ssd " + disk_size + " /cromwell_root/"
 	}
 	output {
-		Array[File] splited_vcf = glob("*.vcf")
+		Array[File] splited_vcf = glob("*.splited.vcf")
 		Array[File] family_vcf = glob("*.family.vcf")
 		File LCL5_1 = "${project}.LCL5_1.vcf"
 		File LCL5_2 = "${project}.LCL5_2.vcf"