SVs impacting genes - mestato/EPP622_2026 GitHub Wiki

Prep

Prepare reference genome. We just need to get a tsv file of contig lengths from the fai file. I have already done this in /lustre/isaac24/proj/UTK0505/module3_data

cut -f1,2 AppSpring_Hap1_Chr9.fasta.fai > AppSpring_Hap1_Chr9.contigLens.tsv

Next prepare the gff3 of genes by sorting, compressing and indexing with tabix. I have already done this in /lustre/isaac24/proj/UTK0505/module3_data.

conda activate svs
module load htslib
(
  grep '^#' AppSpring_Hap1_Chr9.genes.gff3
  grep -v '^#' AppSpring_Hap1_Chr9.genes.gff3 | LC_ALL=C sort -t $'\t' -k1,1 -k4,4n -k5,5n
) > AppSpring_Hap1_Chr9.genes.sorted.gff3
conda install tabix 
bgzip -f AppSpring_Hap1_Chr9.genes.sorted.gff3
tabix -f -p gff AppSpring_Hap1_Chr9.genes.sorted.gff3.gz

run Ensembl VEP

Install ensembl-vep with conda

conda create -n vep ensembl-vep perl-dbi bcftools tabix
conda activate vep

Set up

mkdir 04_vep

Now lets get only our 'PASS' variants from the cuteSV lab and index them .

bcftools view -f PASS ../03_cutesv/CherBrave3.vcf -O z -o CherBrave3.PASS.vcf.gz
tabix -f -p vcf CherBrave3.PASS.vcf.gz

Now we can run VEP. It runs very fast (<1 min) so I'm not even going to bother with slurm.

vep \
    --input_file CherBrave3.PASS.vcf.gz \
    --output_file CherBrave3.PASS.VEP.vcf.gz \
    --vcf --compress_output bgzip --force_overwrite \
    --fasta /lustre/isaac24/proj/UTK0505/module3_data/AppSpring_Hap1_Chr9.fasta \
    --gff /lustre/isaac24/proj/UTK0505/module3_data/AppSpring_Hap1_Chr9.genes.sorted.gff3.gz \
    --distance 1000 

Summary of results can be found in the html file.

Some genes are deleted, lets look at a deletion impacting a gene. Note the new information "transcript_ablation|HIGH":

AppSpring_1.0.0_Hap1_Chr09	94315338	cuteSV.DEL.15332	ATTAATTAAAACCACATG...ATT	A	209.9	PASS	PRECISE;SVTYPE=DEL;SVLEN=-19779;END=94335117;CIPOS=-178,178;CILEN=-8,8;RE=22;AF=1;STRAND=+-;CSQ=-|transcript_ablation|HIGH||g21035|Transcript|g21035.t1|protein_coding|1/1|||||||||||-1||||AppSpring_Hap1_Chr9.genes.sorted.gff3.gz|	GT:DR:DV:PL:GQ	1/1:0:22:210,56,0:56

Get the local bam file

samtools view  -b -h -o CherBrave3.Chr9_94Mb.bam CherBrave3.bam AppSpring_1.0.0_Hap1_Chr09:94265338-94365338
samtools index CherBrave3.Chr9_94Mb.bam

Lets look at a stop gained variant.

CherBrave3.PASS.VEP.vcf:AppSpring_1.0.0_Hap1_Chr09	61736498	cuteSV.INS.9464	C	CTCTTGGTGGGTCTCATTCTATACAACGTACCTCGGCCTCTAAAGCAGGAATTTTAAATATCCAAGGTTGCAATTTTGCAGTGTAAGAGATATATACATCAAATACTTGTAATTCTCTCAATGTAAA	267.2	PASS	PRECISE;SVTYPE=INS;SVLEN=126;END=61736498;CIPOS=-11,11;CILEN=0,0;RE=28;AF=1;CSQ=TCTTGGTGGGTCTCATTCTATACAACGTACCTCGGCCTCTAAAGCAGGAATTTTAAATATCCAAGGTTGCAATTTTGCAGTGTAAGAGATATATACATCAAATACTTGTAATTCTCTCAATGTAAA|stop_gained&inframe_insertion|HIGH||g19730|Transcript|g19730.t1|protein_coding|10/10||||1228-1229|1228-1229|410|L/LLVGLILYNVPRPLKQEF*ISKVAILQCKRYIHQILVILSM*I|ctc/cTCTTGGTGGGTCTCATTCTATACAACGTACCTCGGCCTCTAAAGCAGGAATTTTAAATATCCAAGGTTGCAATTTTGCAGTGTAAGAGATATATACATCAAATACTTGTAATTCTCTCAATGTAAAtc|||1||||AppSpring_Hap1_Chr9.genes.sorted.gff3.gz|	GT:DR:DV:PL:GQ	1/1:0:28:267,71,0:71

Get the local bam file

samtools view  -b -h -o CherBrave3.Chr9_62Mb.bam CherBrave3.bam AppSpring_1.0.0_Hap1_Chr09:61700498-61766498
samtools index CherBrave3.Chr9_62Mb.bam