Merge pull request #7 from CEGRcode/validation

Validation of DelID using YKOC data

Author: owlang

paper/.gitignore

@@ -15,6 +15,10 @@ SyntheticStrain/logs/*.err-*
 SyntheticStrain/logs/*.out-*
 YKOC-wgs/logs/*.err-*
 YKOC-wgs/logs/*.out-*
+YKOC-wgs/results/FASTQ
+YKOC-wgs/results/BAM
+YKOC-wgs/results/ID
+YKOC-wgs/results/BedGraphs
 ENCODEdata-eGFP/logs/*.out-*
 ENCODEdata-eGFP/logs/*.err-*
 ENCODEdata-eGFP/results/FASTQ

paper/YKOC-wgs/210328_Puddu_2019_STable6_verifiedKO.txt

@@ -1,11 +1,15 @@
 # Run DeletionID on WGS data of the YKOC
 
 # Metadata with strain information, EBI accessions, and gene name maps.
-The `201213_Puddu_2019_STable1_del2srs.txt` file was pulled from Supplementary 
+The `201213_Puddu_2019_STable1_del2ers.txt` file was pulled from Supplementary 
 Table 1 from Puddu et al 2019. Maps the EBI sample accessions (ERSXXXXXXX) to the 
 expected knockout strain and a sample ID used by the paper (SDXXXX). This table was 
 downloaded December 13, 2020.
 
+The `210328_Puddu_2019_STable6_verifiedKO.txt` file was pulled from Supplementary 
+Table 6 from Puddu et al 2019. This file contains the list of verified KO samples 
+keyed on their SD id and including associated scores with notes.
+
 The `210403_PRJEB27160_accessions.txt` metadata file was pulled to map EBI sample 
 accessions (ERSXXXXXXX) to the run accessions (ERRXXXXXXX). The file also contains 
 other sequencing information including ftp locations for FASTQ files, MD5 checksums, 
@@ -14,6 +18,14 @@ April 3, 2021 to download this metadata.
 
 https://www.ebi.ac.uk/ena/browser/view/PRJEB27160
 
+Note: Four ERS sample accessions referenced in the paper are missing from the EBI 
+metadata download:
+
+ERS1041599
+ERS969339
+ERS969700
+ERS969701
+
 # Download & Align YKOC data
 Use the EBI accessions to download the data to the `results/FASTQ` directory. The 
 `job/00_download_data.pbs` PBS script uses the `wget` command but for the paper, these

paper/YKOC-wgs/README

@@ -0,0 +1,6 @@
+FLANKING	200
+LOCI	APE3	VAC17	PIR3
+SAMPLES	ERS838258	ERS903113	ERS969760	ERS1076728
+S_LABEL	APE3	VAC17	PIR3	WT-1
+S_COLOR	forestgreen	firebrick	salmon	gray
+S_MAX_Y	200	150	10	150

paper/YKOC-wgs/fig3e_config.txt

@@ -0,0 +1,6 @@
+FLANKING	200
+LOCI	SWT1	PNS1
+SAMPLES	ERS1076684	ERS1076685	ERS1029415	ERS1029416	ERS1076728
+S_LABEL	del_SWT1_rep1	del_SWT1_rep2	del_PNS1_rep1	del_PNS1_rep2	WT-1
+S_COLOR	royalblue	cornflowerblue	goldenrod	gold	gray
+S_MAX_Y	150	150	150	150	150

paper/YKOC-wgs/fig3f_config.txt

@@ -0,0 +1,66 @@
+#!/bin/bash
+#PBS -l nodes=1:ppn=4
+#PBS -l pmem=16gb
+#PBS -l walltime=01:00:00
+#PBS -A open
+#PBS -o logs/align-picard.log.out
+#PBS -e logs/align-picard.log.err
+#PBS -t 1-9010
+
+module load gcc/8.3.1
+module load samtools/1.5
+module load bwa/0.7.15
+module load picard/2.20.8 
+
+WRK=/path/to/GenoPipe/paper/YKOC-wgs
+cd $WRK
+
+[ -d logs ] || mkdir logs
+[ -d results/BAM ] || mkdir results/BAM
+
+GENOME=../input/sacCer3.fa
+METADATA=210403_PRJEB27160_accessions.txt
+INDEX=$(($PBS_ARRAYID+1))
+
+INFO=`sed "${INDEX}q;d" $METADATA`
+ERR=`echo $INFO | awk '{print $4}'`
+ERS=`echo $INFO | awk '{print $2}'`
+FQ1=results/FASTQ/$ERR/$ERR\_1.fastq.gz
+FQ2=results/FASTQ/$ERR/$ERR\_2.fastq.gz
+BAM=$WRK/results/BAM/$ERS
+#echo $INFO
+
+start=`date +%s`
+# Align with Pugh Lab standard alignment pipeline parameters
+# -T INT  Don't output alignments with score lower than INT. This option affects output and occsaionally SAM flag 2.
+# -h INT  If query has not more than INT hits with score higher than 880% of the best hit, output them all in the XA tag.
+# -M      Mark shorter split hits as secondary (for Picard compatibility).
+echo "(${PBS_ARRAYID}) Aligning $ERR fastq files..."
+bwa mem -v 1 -T '30' -h '5' -t 4 -M $GENOME $FQ1 $FQ2 > $BAM.unsorted.bam
+# Sort
+echo "(${PBS_ARRAYID}) Sorting $ERR ..."
+samtools sort $BAM.unsorted.bam > $BAM.unmarked.bam
+# Mark duplicates
+echo "(${PBS_ARRAYID}) Marking duplicates $ERR ..."
+picard MarkDuplicates \
+	INPUT=$BAM.unmarked.bam \
+	OUTPUT=$BAM.marked.bam \
+	METRICS_FILE=$BAM.metrics.txt \
+	REMOVE_DUPLICATES='false' ASSUME_SORTED='true' \
+	DUPLICATE_SCORING_STRATEGY='SUM_OF_BASE_QUALITIES' \
+	#READ_NAME_REGEX='[a-zA-Z0-9]+:[0-9]:([0-9]+):([0-9]+):([0-9]+).*.' \
+	#OPTICAL_DUPLICATE_PIXEL_DISTANCE='100' \
+	#VALIDATION_STRINGENCY='LENIENT' #VERBOSITY=ERROR
+# Filter duplicates
+echo "(${PBS_ARRAYID}) Filter duplicates $ERR ..."
+samtools view  -h -b -f 0x1 -F 0x404 \
+	-o $BAM.bam \
+	$BAM.marked.bam
+# Index
+echo "(${PBS_ARRAYID}) Index $ERR ..."
+samtools index $BAM.bam
+end=`date +%s`
+runtime=$((end-start))
+echo "...alignment and indexing for ($PBS_ARRAYID) $ERR.fq > $ERS.bam finished in ${runtime}"
+# Clean-up
+rm $BAM.unsorted.bam $BAM.unmarked.bam $BAM.marked.bam

paper/YKOC-wgs/job/01_align_fastq.pbs

@@ -0,0 +1,59 @@
+#!/bin/bash
+#PBS -l nodes=1:ppn=4
+#PBS -l pmem=16gb
+#PBS -l walltime=00:10:00
+#PBS -A open
+#PBS -o logs/did-picard.log.out
+#PBS -e logs/did-picard.log.err
+#PBS -t 1-9010
+
+module load gcc/8.3.1
+module load anaconda3
+source activate genopipe
+
+# Python3.9.0 env with pysam
+
+WRK=/path/to/GenoPipe/paper/YKOC-wgs
+cd $WRK
+
+METADATA=210403_PRJEB27160_accessions.txt
+INDEX=$(($PBS_ARRAYID+1))
+
+INFO=`sed "${INDEX}q;d" $METADATA`
+ERS=`echo $INFO | awk '{print $2}'`
+BAM=$WRK/results/BAM/$ERS
+ID=$WRK/results/ID
+TEMP=$WRK/temp$PBS_ARRAYID
+#echo $INFO
+
+[ -d logs ] || mkdir logs
+[ -d $ID ] || mkdir $ID
+
+# Check if ID already generated...
+#if [[ -f $ID/$ENCFF\_R1-ID.tab ]]; then
+#	echo "ID already generated ($ENCFF). Exiting.."
+#	exit
+#fi
+
+[ -d $TEMP ] || mkdir $TEMP
+[ -d $ID ] || mkdir -p $ID
+
+# Set-up Temp directory
+cd $TEMP
+ln -s $BAM.bam
+ln -s $BAM.bam.bai
+
+DATABASE=$WRK/../db/sacCer3_Del
+GENOPIPE=$WRK/../../
+
+## Execute Mass EpitopeID and record time
+cd $WRK
+cd $GENOPIPE/DeletionID
+echo "(2) Begin executing EpitopeID for $R1..."
+start=`date +%s`
+bash identify-Deletion.sh -i $TEMP -o $ID -d $DATABASE
+end=`date +%s`
+runtime=$((end-start))
+echo "...single DeletionID for ($PBS_ARRAYID) $R1 finished in ${runtime}"
+
+rm -r $TEMP

paper/YKOC-wgs/job/02_indexed_runDID.pbs

@@ -0,0 +1,39 @@
+#!/bin/bash
+#PBS -l nodes=1:ppn=4
+#PBS -l pmem=16gb
+#PBS -l walltime=00:10:00
+#PBS -A open
+#PBS -o tally-results.out
+#PBS -e tally-results.err
+
+module load anaconda3
+source activate genopipe
+#PyLib argparse
+#PyLib matplotlib
+#PyLib matplotlib-venn
+
+# FIRST CHANGE PATH TO EXECUTE
+WRK=/path/to/GenoPipe/paper/YKOC-wgs
+cd $WRK
+
+#wget http://sgd-archive.yeastgenome.org/curation/chromosomal_feature/SGD_features.tab
+wget http://sgd-archive.yeastgenome.org/curation/chromosomal_feature/deleted_merged_features.tab
+wget http://sgd-archive.yeastgenome.org/curation/chromosomal_feature/saccharomyces_cerevisiae.gff.gz
+gzip -d saccharomyces_cerevisiae.gff.gz
+
+# Tally results to generate TableS2-ykoc-output-summary
+python scripts/analyze_ykoc_results.py -m 201213_Puddu_2019_STable1_del2ers.txt -i results/ID \
+	-e 210403_PRJEB27160_accessions.txt -t 210328_Puddu_2019_STable6_verifiedKO.txt \
+	-d deleted_merged_features.tab -n saccharomyces_cerevisiae.gff \
+	> results/ykoc_output_summary.txt \
+	2> results/ykoc_output_summary.err
+
+grep 'PASS' results/ykoc_output_summary.txt | grep $'\tVerified' > results/ConfidentPass.txt
+grep 'PASS' results/ykoc_output_summary.txt | grep 'NotVerified' > results/TentativePass.txt
+grep 'FAIL' results/ykoc_output_summary.txt | grep $'\tVerified' > results/TentativeFails.txt
+grep 'FAIL' results/ykoc_output_summary.txt | grep 'NotVerified' > results/TrueFails.txt
+
+python scripts/make_venn_diagram.py -s results/ykoc_output_summary.txt -t Figure3C
+
+#Clean-up
+rm saccharomyces_cerevisiae.gff* deleted_merged_features.tab

paper/YKOC-wgs/job/03_tally_results.pbs

@@ -0,0 +1,46 @@
+#!/bin/bash
+#PBS -l nodes=1:ppn=4
+#PBS -l pmem=16gb
+#PBS -l walltime=03:00:00
+#PBS -A open
+#PBS -o logs/pileup.samples.log.out
+#PBS -e logs/pileup.samples.log.err
+#PBS -t 1-9010
+
+module load gcc/8.3.1
+module load samtools/1.5
+module load bedtools/2.27.1
+module load anaconda3
+source activate genopipe
+
+WRK=/path/to/GenoPipe/paper/YKOC-wgs
+cd $WRK
+
+[ -d logs ] || mkdir logs
+[ -d results/BedGraphs ] || mkdir results/BedGraphs
+
+METADATA=results/ykoc_output_summary.txt
+INDEX=$(($PBS_ARRAYID+1))
+
+INFO=`sed "${INDEX}q;d" $METADATA`
+ERS=`cut -f10 $METADATA |sed "${INDEX}q;d"`
+BAM=results/BAM/$ERS.bam
+BG=results/BedGraphs/$ERS
+#echo $INFO
+
+# Pileup to BedGraph
+echo "($PBS_ARRAYID) Make raw bedgraph pileups for $ERS"
+bedtools genomecov -bga -scale 1.0 -ibam $BAM \
+	| sort -k1,1 -k2,2n \
+	> $BG.raw.bedgraph
+
+# Normalize BedGraph
+echo "($PBS_ARRAYID) Normalize bedgraph for $ERS"
+NORMALIZE=scripts/normalize_BedGraph.py
+python $NORMALIZE -i $BG.raw.bedgraph > $BG.bedgraph
+
+# Make CDT row centered on expected KO
+echo "($PBS_ARRAYID) Generate CDT slice for $ERS"
+ROWCDT=scripts/make_cdt_row.py
+echo $INFO
+python $ROWCDT -s <(head -n $INDEX $METADATA |tail -n 1) -i results/BedGraphs/ -g saccharomyces_cerevisiae.gff

paper/YKOC-wgs/job/04_pileup_samples.pbs

@@ -0,0 +1,71 @@
+#!/bin/bash
+#PBS -l nodes=1:ppn=4
+#PBS -l pmem=16gb
+#PBS -l walltime=05:00:00
+#PBS -A open
+#PBS -o logs/make.figs.log.out
+#PBS -e logs/make.figs.log.err
+
+#1-9010
+#bfp2_h_g_sc_default
+#bfp2_j_g_bc_default
+
+module load gcc/8.3.1
+module load samtools/1.5
+module load bedtools/2.27.1
+module load anaconda3
+source activate genopipe
+
+WRK=/path/to/GenoPipe/paper/YKOC-wgs
+WRK=/storage/home/owl5022/scratch/GenotypingProject/GenoPipe/paper/YKOC-wgs
+cd $WRK
+
+[ -d logs ] || mkdir logs
+[ -d results/figs ] || mkdir results/figs
+
+for PARTITION in "TentativePass" "TrueFails" "TentativeFails" "ConfidentPass";
+do
+	BASE=results/$PARTITION
+	rm $BASE.temp
+	
+	METADATA=results/$PARTITION.txt
+	LNCT=`wc -l $METADATA | awk '{print $1}'`
+	for ((INDEX=1; INDEX<=$LNCT; INDEX++));
+	do
+		#echo $INDEX
+		ERS=`cut -d $'\t' -f10 $METADATA | sed "${INDEX}q;d"`
+		ORF=`cut -d $'\t' -f7 $METADATA | sed "${INDEX}q;d"`
+		echo "($INDEX) Parsed out ERS=$ERS and ORF=$ORF"
+		ROW=results/BedGraphs/$ERS\_$ORF\_6000bp.cdt
+		cat $ROW >> $BASE.temp
+	done
+	
+	sort -nk3,3 $BASE.temp > $BASE.cdt
+	
+	# Make Heatmap PNG
+	##ScriptManager Two-color heatmap
+	# Black, .95 threshold
+	# $BASE.cdt > $BASE.png
+
+	# Make Trace CDT
+	TRACE=scripts/make_trace.py
+	python $TRACE -c $BASE.cdt -g saccharomyces_cerevisiae.gff -w 6000 > $BASE\_trace.cdt
+
+	rm $BASE.temp	
+done
+
+# Make trace (x4 lines)
+##Script to make trace cdt
+##ScriptManager Three-color heatmap
+# diff colors, set alpha channel, threshold 0<0<.95
+
+SUBPLOTS=scripts/make_subplots.py
+
+# Figure 3E SVG image
+python $SUBPLOTS -t Figure3E -c fig3e_config.txt -g saccharomyces_cerevisiae.gff
+
+# Figure 3F SVG image
+python $SUBPLOTS -t Figure3F -c fig3f_config.txt -g saccharomyces_cerevisiae.gff
+
+
+