Skip to content

Remove b37 dna #95

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Istaisa wants to merge 3 commits into
base: development
Choose a base branch
from
Open

Remove b37 dna #95

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
27971d5
Remove b37 TSO500 DNA
Istaisa Apr 13, 2026
03ac1e6
update readme
Istaisa Apr 13, 2026
5db8a7d
remove unnecessary conda environment and activation
Istaisa Apr 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
278 changes: 10 additions & 268 deletions 2_TSO500.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,8 @@
#SBATCH --partition=high
#SBATCH --cpus-per-task=24

# Description: Run Illumina TSO500 app for each sample then run postprocessing steps - FastQC,
# GATK depth of coverage, coverage calculator, bed2hgvs, gather QC metrics. Kick
# off script 3 when all samples completed
# Description: Run Illumina TSO500 app for each sample then run postprocessing steps for RNA only.
# Kick off script 3 when all samples completed
# Use: from /Output/results/<run_id>/TSO500/ directory, for each sample run:
# sbatch --export=raw_data=/data/raw/novaseq/<run_id>,sample_id=<sample_id> 2_TSO500.sh
# Version: 1.0.15
Expand Down Expand Up @@ -40,15 +39,6 @@ module load anaconda
# catch fails early and terminate
set -euo pipefail

# define pipeline variables
minimum_coverage="270 135"
coverage_bed_files_path="$pipeline_dir"/hotspot_coverage
vendor_capture_bed="$pipeline_dir"/vendorCaptureBed_100pad_updated.bed
preferred_transcripts="$pipeline_dir"/preferred_transcripts.txt
worksheet=$(grep "$sample_id" SampleSheet_updated.csv | cut -d, -f2)
dna_or_rna=$(grep ${sample_id}, SampleSheet_updated.csv | cut -d, -f9)


##############################################################################################
# Illumina app
##############################################################################################
Expand Down Expand Up @@ -99,219 +89,10 @@ for fastqPair in $(find $fastq_path -name *.fastq.gz -type f -printf "%f\n" | cu
done


##############################################################################################
# DNA only steps
##############################################################################################

if [ "$dna_or_rna" = "DNA" ]; then


#-------------------------------------------------------------------------------------
# Call variants outside of app ROI
#-------------------------------------------------------------------------------------

bash "$pipeline_scripts"/call_extra_padding_variants.sh "$sample_id" "$pipeline_version"


#-------------------------------------------------------------------------------------
# Run depth of coverage with limited bed (whole 500 takes a long time)
#-------------------------------------------------------------------------------------

# set and make filepaths for depth of coverage
bam_path="$output_path"/Logs_Intermediates/StitchedRealigned/"$sample_id"/
depth_path="$output_path"/depth_of_coverage
mkdir -p $depth_path

# reheader the bams to local area
java -jar /Apps/wren/picard/2.21.6/bin/picard.jar AddOrReplaceReadGroups \
I="$bam_path"/"$sample_id".bam \
O="$bam_path"/"$sample_id"_add_rg.bam \
RGID=4 \
RGLB=lib1 \
RGPL=ILLUMINA \
RGPU=unit1 \
RGSM=20

# index new bam
samtools index "$bam_path"/"$sample_id"_add_rg.bam "$bam_path"/"$sample_id"_add_rg.bam.bai

# run depth of coverage
gatk DepthOfCoverage \
-I "$bam_path"/"$sample_id"_add_rg.bam \
-L "$vendor_capture_bed" \
-R "$app_dir"/resources/genomes/hg19_hardPAR/genome.fa \
-O "$depth_path"/"$sample_id"_depth_of_coverage

# change to tab delimited and remove colon from column 1
sed 's/:/\t/g' "$depth_path"/"$sample_id"_depth_of_coverage \
| sed 's/,/\t/g' | grep -v 'Locus' \
| sort -k1,1 -k2,2n | bgzip \
> "$depth_path"/"$sample_id"_depth_of_coverage.gz

# tabix index depth of coverage file
tabix \
-b 2 \
-e 2 \
-s 1 \
"$depth_path"/"$sample_id"_depth_of_coverage.gz

# deactivate env
set +u
conda deactivate
set -u


#-------------------------------------------------------------------------------------
# Run coverage calculator at 135X and 270X depth cutoffs
#-------------------------------------------------------------------------------------

# repeat for each coverage value
for min_coverage in $minimum_coverage; do

# activate coverage calculator conda env
set +u
conda activate CoverageCalculatorPy
set -u

# set output directory for coverage files
hscov_outdir=hotspot_coverage_"$min_coverage"x

# run coverage calculator on each bed file
for bed_file in "$coverage_bed_files_path"/*.bed; do

name=$(echo $(basename $bed_file) | cut -d"." -f1)

python /data/diagnostics/apps/CoverageCalculatorPy/CoverageCalculatorPy-v1.1.0/CoverageCalculatorPy.py \
-B "$coverage_bed_files_path"/"$name".bed \
-D "$depth_path"/"$sample_id"_depth_of_coverage.gz \
--depth "$min_coverage" \
--padding 0 \
--groupfile "$coverage_bed_files_path"/"$name".groups \
--outname "$sample_id"_"$name" \
--outdir "$depth_path"/"$hscov_outdir"/

# remove header from gaps file
if [[ $(wc -l < "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".gaps) -eq 1 ]]; then
# no gaps
touch "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".nohead.gaps

else
# gaps
grep -v '^#' "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".gaps > "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".nohead.gaps

fi

# remove chr from bed file so bed2hgvs works
cat "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".nohead.gaps | sed 's/^chr//' > "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".nohead_nochr.gaps

# remove intermediate files
rm "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".gaps
rm "$depth_path"/"$hscov_outdir"/"$sample_id"_"$name".nohead.gaps

done


#-------------------------------------------------------------------------------------
# Run bed2hgvs to add hgvs nomenclature to gaps
#-------------------------------------------------------------------------------------

# activate bed2hgvs conda env
set +u
conda deactivate
conda activate bed2hgvs
set -u

# run on each bed file
for gaps_file in "$depth_path"/"$hscov_outdir"/*.nohead_nochr.gaps; do

name=$(echo $(basename $gaps_file) | cut -d"." -f1)
echo $name

Rscript /data/diagnostics/apps/bed2hgvs/bed2hgvs-v0.3.0/bed2hgvs.R \
--bedfile $gaps_file \
--outname "$name".gaps \
--outdir "$depth_path"/"$hscov_outdir" \
--preferred_tx $preferred_transcripts

# remove intermediate file
rm "$depth_path"/"$hscov_outdir"/"$name".nohead_nochr.gaps
done

# combine all total coverage files
if [ -f "$depth_path"/"$hscov_outdir"/"$sample_id"_coverage.txt ]; then rm "$depth_path"/"$hscov_outdir"/"$sample_id"_coverage.txt; fi
cat "$depth_path"/"$hscov_outdir"/*.totalCoverage | grep "FEATURE" | head -n 1 >> "$depth_path"/"$hscov_outdir"/"$sample_id"_coverage.txt
cat "$depth_path"/"$hscov_outdir"/*.totalCoverage | grep -v "FEATURE" | grep -vP "combined_\\S+_GENE" >> "$depth_path"/"$hscov_outdir"/"$sample_id"_coverage.txt

# deactivate env
set +u
conda deactivate
set -u


#-------------------------------------------------------------------------------------
# Cosmic gaps
#-------------------------------------------------------------------------------------

# activate conda env
set +u
conda activate TSO500_post_processing
set -u

cosmic_tool_path=/data/diagnostics/apps/cosmic_gaps/cosmic_gaps-master

# parse referral - must be in DNA loop
referral=$(grep "$sample_id" samples_correct_order_"$worksheet"_DNA.csv | cut -d, -f4)
gaps_file="$depth_path"/"$hscov_outdir"/"$sample_id"_"$referral"_hotspots.gaps

# hotspot gaps file may be missing for some referrals
if [[ -f $gaps_file ]]
then

# only run bedtools intersect for certain referral types
if [ $referral = "Melanoma" ] || [ $referral = "Lung" ] || [ $referral = "Colorectal" ] || [ $referral = "GIST" ] || [ $referral = "breast" ]
then
dos2unix $gaps_file

# find the overlap between the hotspots file and the referral file from cosmic
bedtools intersect \
-loj \
-F 1 \
-a $gaps_file \
-b "$cosmic_tool_path"/cosmic_bedfiles/"$referral".bed \
-wao \
> "$depth_path"/"$hscov_outdir"/"$sample_id"_"$referral"_intersect.txt

fi

# filter the output
python "$cosmic_tool_path"/filter_table.py \
--sampleId $sample_id \
--referral $referral \
--gaps_path "$depth_path"/"$hscov_outdir"/ \
--bedfile_path "$cosmic_tool_path"/cosmic_bedfiles/

fi

# deactivate env
set +u
conda deactivate
set -u

done

fi



##############################################################################################
# Gather QC metrics for sample
##############################################################################################

# activate conda env
set +u
conda activate TSO500_post_processing
set -u

# function to check FASTQC output
count_qc_fails() {
#count how many core FASTQC tests failed
Expand All @@ -334,54 +115,15 @@ done
# pull out metrics from the Illumina app MetricsOutput.tsv
completed_all_steps=$(grep COMPLETED_ALL_STEPS analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f2)

# DNA only metrics
if [ "$dna_or_rna" = "DNA" ]; then

contamination_score=$(grep CONTAMINATION_SCORE analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f4)
contamination_p_value=$(grep CONTAMINATION_P_VALUE analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f4)
total_pf_reads=$(grep TOTAL_PF_READS analysis/"$sample_id"/Results/MetricsOutput.tsv | head -n1 | cut -f4)
median_insert_size=$(grep MEDIAN_INSERT_SIZE analysis/"$sample_id"/Results/MetricsOutput.tsv | head -n1 | cut -f4)
median_exon_coverage=$(grep MEDIAN_EXON_COVERAGE analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f4)
pct_exon_50x=$(grep PCT_EXON_50X analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f4)

contamination_pass_fail=PASS
if [ $contamination_score -gt 3106 ]; then

# use python to handle float
p_value_cutoff=$(python -c "print(float("$contamination_p_value") > 0.049)")
if [ $p_value_cutoff = "True" ]; then
contamination_pass_fail=FAIL
fi
fi


#run samtools to get the number of reads in the bam file
if [[ -f ./analysis/"$sample_id"/Logs_Intermediates/StitchedRealigned/"$sample_id"/"$sample_id".bam ]]; then
reads=$( samtools view -c ./analysis/"$sample_id"/Logs_Intermediates/StitchedRealigned/"$sample_id"/"$sample_id".bam )
else
reads="NA"
fi

# add to sample QC file
echo -e "Sample\tFastQC\tcompleted_all_steps\tcontamination_pass_fail\tcontamination_score\tcontamination_p_value\ttotal_pf_reads\tmedian_insert_size\tmedian_exon_coverage\tpct_exon_50x\tAligned_reads" > "$output_path"/"$sample_id"_"$dna_or_rna"_QC.txt
echo -e "$sample_id\t$fastqc_status\t$completed_all_steps\t$contamination_pass_fail\t$contamination_score\t$contamination_p_value\t$total_pf_reads\t$median_insert_size\t$median_exon_coverage\t$pct_exon_50x\t$reads" >> "$output_path"/"$sample_id"_"$dna_or_rna"_QC.txt

fi


# RNA only metrics
if [ "$dna_or_rna" = "RNA" ]; then

median_cv_gene_500x=$(grep "MEDIAN_CV_GENE_500X" analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f4)
total_on_target_reads=$(grep "TOTAL_ON_TARGET_READS" analysis/"$sample_id"/Results/MetricsOutput.tsv | tail -n1 | cut -f4)
median_insert_size=$(grep "MEDIAN_INSERT_SIZE" analysis/"$sample_id"/Results/MetricsOutput.tsv | tail -n1 | cut -f4)
total_pf_reads=$(grep "TOTAL_PF_READS" analysis/"$sample_id"/Results/MetricsOutput.tsv | tail -n1 | cut -f4)

# add to sample QC file
echo -e "Sample\tFastQC\tcompleted_all_steps\tmedian_cv_gene_500x\ttotal_on_target_reads\tmedian_insert_size\ttotal_pf_reads" > "$output_path"/"$sample_id"_"$dna_or_rna"_QC.txt
echo -e "$sample_id\t$fastqc_status\t$completed_all_steps\t$median_cv_gene_500x\t$total_on_target_reads\t$median_insert_size\t$total_pf_reads" >> "$output_path"/"$sample_id"_"$dna_or_rna"_QC.txt

fi
median_cv_gene_500x=$(grep "MEDIAN_CV_GENE_500X" analysis/"$sample_id"/Results/MetricsOutput.tsv | cut -f4)
total_on_target_reads=$(grep "TOTAL_ON_TARGET_READS" analysis/"$sample_id"/Results/MetricsOutput.tsv | tail -n1 | cut -f4)
median_insert_size=$(grep "MEDIAN_INSERT_SIZE" analysis/"$sample_id"/Results/MetricsOutput.tsv | tail -n1 | cut -f4)
total_pf_reads=$(grep "TOTAL_PF_READS" analysis/"$sample_id"/Results/MetricsOutput.tsv | tail -n1 | cut -f4)

# add to sample QC file
echo -e "Sample\tFastQC\tcompleted_all_steps\tmedian_cv_gene_500x\ttotal_on_target_reads\tmedian_insert_size\ttotal_pf_reads" > "$output_path"/"$sample_id"_RNA_QC.txt
echo -e "$sample_id\t$fastqc_status\t$completed_all_steps\t$median_cv_gene_500x\t$total_on_target_reads\t$median_insert_size\t$total_pf_reads" >> "$output_path"/"$sample_id"_RNA_QC.txt


##############################################################################################
Expand Down
44 changes: 11 additions & 33 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
# TSO500_post_processing

This pipeline
* Demultiplexes with the Illumina app
* Runs the RNA samples through the Illumina app
* Submits the DNA samples to run through the [somatic_enrichment_nextflow](https://github.com/AWGL/somatic_enrichment_nextflow) pipeline

## Documentation

Expand All @@ -16,30 +20,23 @@ From this folder run the command:

The raw data directory must contain the SampleSheet.csv.

## To re-run DNA samples with a new referral
## To run a DNA sample

The script rerun_coverage.sh allows the coverage json file to be produced for a DNA sample with a new referral in all scenarios (analysis/ folder present or deleted). This can be run as follows, as transfer, from any location:

`sbatch rerun_coverage.sh <run_id> <sample_id> <worksheet> <new referral>`

The new coverage file will then be present in Gathered_Results/Database/ alongside a new samples list file for upload to the database.
To run a DNA sample through the Illumina app, please run an old release (v1.1.1)

## Duty scientist responsibilites
The duty scientist is responsible for the following tasks:
* Creating the TSO500 samplesheet- see below for requirements
* Signing off the run in autoqc database - for more information refer to the AutoQC sop
* Importing the data into the data into the somatic variant database

* Signing off the run in autoqc database and GLIMS (NOTE: the RNA must now be signed off with the DNA)
* Signing off the run in SVD

## Samplesheet requirements
* The samplesheet must contain the samples in the correct order for the RNA contamination results to be valid
* Every NTC must be named NTC-worksheetid



## Unit tests

Unit tests have been created against the following scripts: `tsv2db.py`, `coverage2json.py`, `fusion_check_with_ntc.py`.
Unit tests have been created against the following scripts: `fusion_check_with_ntc.py`.

To run all unit tests:
- copy these scripts into the `tests/` folder (relative imports not currently working, will be fixed in future version)
Expand All @@ -50,34 +47,15 @@ To run all unit tests:
To run tests on a specific script, follow the steps above but run `python -m unittest <test_script_name>`


## Adding a new DNA panel

To add a new panel, the following needs to be changed:

**Samplesheet generator:**
- Add the referral reason to the samplesheet generator (see SOP in Qpulse) and make sure it matches the filename of the new bed files (case sensitive)

**Pipeline:**
- Generate bed files for the new panel
- The `hotspot_variants/*bed` bed file and `hotspot_coverage/*combined.bed` files are required
- The `hotspot_coverage/*hotspots.bed` and `hotspot_coverage/*genescreen.bed` files are optional
- Filenames should be all lowercase
- Make sure that all regions in the new panel are covered in `vendorCaptureBed_100pad_updated.bed` (the bed file file used to generate the depth of coverage file)
- Make sure that any flanking regions are added to the `TSO_extra_padding_chr.interval_list` file - Illumina bed file only goes +/- 2bp so this file contains the extra 3bp to make it +/- 5bp

**Somatic variant database:**
- Make a new panel object in the somatic variant database that matches the filename of the new bed files (case sensitive)
- Move the new variants bed file into the `roi/variant_calling` folder in the somatic database

## Adding a new RNA panel

To add a new panel, the following needs to be changed:

**Samplesheet generator:**
- Add the referral reason to the samplesheet generator (see SOP in Qpulse)
- Add the referral to the SampleSheet Generator, including the mapped test directory code

**Pipeline:**
- Add a new file to RNA_referrals named <panel>.txt with the gene names on the panel, one per line

**Somatic variant database:**
- Make a new panel object in the somatic variant database that includes the genes on the panel to filter fusions by
- Make a new Panel object in SVD
Loading