ahmedarslan/doppelganger

Target site duplication assessment from alignment file and genomic location of insertion sites.

doppelganger

Target site duplication assessment from alignment file and genomic location of somatic insertion sites.

Instructions to generate input files:

• There are two ways to generate input files.

  1 - follow individual steps - OR

  2 - run the preprocessor.sh to automate files generation (Recommended)

• Three files are needed:

  1 - bam file (sorted)

  2 - bed file (containing insertion locations)

  3 - fasta file (containing insertion sequences)

1 - follow individual steps

• Step - 1 - Generating (sorted) BAM file: Softwares/data needed: minimap2, samtools, reference

Alignment:
Command: minimap2 -ax map-hifi -MD GRCh38_genomic.fa raw.fastq.gz > HiFiCCS.sam -t 10

  A - change map-hifi to map-ont, if using nanopore reads.
  B - change -t to increase or decrease number of threads. 
  C - change reference genome (GRCh38_genomic.fa) to your requirement.

Convert sam to bam file format:
Command: samtools view -Sb HiFiCCS.sam -@ 10 > HiFiCCS.bam 
  A - change -@ to increase or decrease number of threads.

Sort bam file:
Command: samtools sort -o HiFiCCS.sorted.bam -O bam -@ 10 HiFiCCS.bam 
Index sorted bam file:
Command: samtools index HiFiCCS.sorted.bam

• Step - 2 - Generating bed and fasta from bed file: Softwares needed: Sniffles, bcftools

Step - Insertion variant calling:
Command: sniffles --input HiFiCCS.sorted.bam --vcf HiFiCCS.sorted.vcf --threads 10 --reference GRCh38_genomic.fa --non-germline --minsupport 1

  A - change reference genome (GRCh38_genomic.fa) to your requirement.
  B - change to germline variant prediction by removing option --non-germline
  C - increase number of reads required to predict an insertion by changing option --minsupport [int]
  
Step - vcf to bed/fasta:
Command: bcftools view -i 'GT="alt"' -f PASS -c 1 HiFiCCS.sorted.vcf -o HiFiCCS.filter.vcf
Command: bcftools query -f '%CHROM\t%POS\t%INFO/END\t%INFO/SVTYPE\t%INFO/SVLEN\t%REF\t%ALT\t%STRAND\n' HiFiCCS.filter.vcf > HiFiCCS.filter.bed
Command: grep 'INS' HiFiCCS.filter.bed > HiFiCCS.filter.INS.bed
Command: awk '{print ">"$1":"$2"-"$3"_"$5"\n"$7}' HiFiCCS.filter.INS.bed > HiFiCCS.fa
Command: awk '{print ">"$1":"$2"-"$3"_"$4"_"$5}' HiFiCCS.filter.INS.bed | awk -F'[->_:]' '{print $2"\t"$3"\t"$3+$6}' | sort -k 1,1 -k2,2n > HiFiCCS.bed

    * HiFiCCS.fa, HiFiCCS.bed and HiFiCCS.sorted.bam are the input files for doppelganger. 

Remove extra/tmp files: rm HiFiCCS.sam HiFiCCS.bam HiFiCCS.sorted.vcf HiFiCCS.filter.bed HiFiCCS.filter.INS.bed HiFiCCS.filter.vcf

2 - run the preprocessor.sh to automate files generation

• conda env create --file doppelganger.yaml

• preprocessor.sh [fastq] [barcode.fa] [threads] [out_name] [reference]

  fastq = raw long-reads file [fastq/fq] ; barcode.fa = pacbio barcode [fa]; out_name = name of output file; reference = refernece genome [fa]

run doppelganger

• conda env create --file doppelganger.yaml [run first time only]
• activate doppelganger
• python3 doppelganger.py --bam [bam file] --bed [bed file] --out [output name] --region [region around insertion] --org [human or mouse] --fasta [insertion sequences.fa] --threads [number of threads]