catching up

Author: ACEnglish

realignment/scripts/get_reference.py

@@ -4,10 +4,15 @@
 
 
 fasta = pysam.FastaFile(sys.argv[1])
-chrom, rest = sys.argv[2].split(':')
-start, end = rest.split('-')
-start = int(start)
-end = int(end)
+if ':' in sys.argv[2]:
+    chrom, rest = sys.argv[2].split(':')
+    start, end = rest.split('-')
+    start = int(start)
+    end = int(end)
+else:
+    chrom = sys.argv[2]
+    start = 1
+    end = fasta.get_reference_length(chrom)
 
 oseq = fasta.fetch(chrom, start - 1, end)
 oseq = re.sub("(.{60})", "\\1\n", oseq, 0, re.DOTALL)

realignment/scripts/msa_realign_pipe.sh

@@ -12,11 +12,14 @@ echo "MSA realignment of region" $region
 mkdir -p msa_$region
 
 #bcftools norm -m - -c s -f $ref -r $region $vcf > msa_${region}/variants.vcf
-bcftools view -r $region -o msa_${region}/variants.vcf.gz -O z $vcf
+bcftools view -r $region $vcf \
+    | bcftools +fill-from-fasta /dev/stdin -- -c REF -f $ref \
+    | bgzip > msa_${region}/variants.vcf.gz
+
 tabix msa_${region}/variants.vcf.gz
 
 python $DIR/get_reference.py $ref $region > msa_${region}/haps.fa
-for i in $(zgrep -m1 '#CHROM' $vcf | cut -f10-)
+for i in $(bcftools view -h $vcf | grep -m1 '#CHROM' | cut -f10-)
 do
     samtools faidx $ref $region | bcftools consensus -H1 --sample $i $vcf | python $DIR/fa_rename.py ${i}_1 >> msa_${region}/haps.fa
     samtools faidx $ref $region | bcftools consensus -H2 --sample $i $vcf | python $DIR/fa_rename.py ${i}_2 >> msa_${region}/haps.fa
@@ -25,7 +28,7 @@ python $DIR/remove_redundant.py msa_${region}/haps.fa > msa_${region}/haps_nored
 
 
 #/users/u233287/scratch/misc_software/mafft-linux64/mafft.bat --auto msa_${region}/haps_noredund.txt > msa_${region}/aln_results.txt
-/users/u233287/scratch/misc_software/mafft-linux64/mafft.bat --auto msa_${region}/haps.fa > msa_${region}/aln_results.txt
+/users/u233287/scratch/misc_software/mafft-linux64/mafft.bat --retree 2 --maxiterate 0 msa_${region}/haps.fa > msa_${region}/aln_results.txt
 #/users/u233287/scratch/misc_software/mafft-linux64/mafft.bat --globalpair --maxiterate 1000 msa_${region}/haps_noredund.txt > msa_${region}/aln_results.txt
 
 #./ProGraphMSA+TR.sh -o result_${region}.txt -R haps_noredund_${region}.txt
@@ -42,17 +45,18 @@ bash $DIR/mend_report.sh msa_${region}/variants.vcf.gz >> msa_${region}/report.t
 echo "Realigned" >> msa_${region}/report.txt
 bash $DIR/mend_report.sh msa_${region}/result.vcf.gz >> msa_${region}/report.txt
 
-python $DIR/get_reference.py $ref $region > msa_${region}/haps_final.fa
-for i in $(zgrep -m1 '#CHROM' $vcf | cut -f10-)
-do
-    samtools faidx $ref $region \
-        | bcftools consensus -H1 --sample $i msa_${region}/result.vcf.gz \
-        | python $DIR/fa_rename.py ${i}_1 >> msa_${region}/haps_final.fa
-    samtools faidx $ref $region \
-        | bcftools consensus -H2 --sample $i msa_${region}/result.vcf.gz \
-        | python $DIR/fa_rename.py ${i}_2 >> msa_${region}/haps_final.fa
-done
-
-echo "md5sums" $(sort msa_${region}/haps.fa | md5sum) $(sort msa_${region}/haps_final.fa | md5sum)
+# Turning off validation checking for now
+#python $DIR/get_reference.py $ref $region > msa_${region}/haps_final.fa
+#for i in $(zgrep -m1 '#CHROM' $vcf | cut -f10-)
+#do
+    #samtools faidx $ref $region \
+        #| bcftools consensus -H1 --sample $i msa_${region}/result.vcf.gz \
+        #| python $DIR/fa_rename.py ${i}_1 >> msa_${region}/haps_final.fa
+    #samtools faidx $ref $region \
+        #| bcftools consensus -H2 --sample $i msa_${region}/result.vcf.gz \
+        #| python $DIR/fa_rename.py ${i}_2 >> msa_${region}/haps_final.fa
+#done
+#
+#echo "md5sums" $(sort msa_${region}/haps.fa | md5sum) $(sort msa_${region}/haps_final.fa | md5sum)
 
 

regions/intersection/README.md

@@ -89,3 +89,27 @@ So, 3.81% of the genome which the v0.2 TR regions cover contains
 - 45.2% of all variants by bases effected
 - 75.5% of SVs by count
 - 47.0% of SVs by bases effected
+
+
+Question 3
+==========
+Can we find expansions/contractions of the tr_annotations inside the variants?
+
+The `filtered_variants_to_regions.txt` is now our new version of the tr_regions.bed. We'll use that to repeat the
+'Defining Repeats' steps described in `../README.md`
+
+
+```bash
+samtools faidx -r <(zcat tr_regions.bed.gz | awk '{print $1 ":" $2 "-" $3}')
+~/scratch/insertion_ref/msru/data/reference/grch38/GRCh38_1kg_mainchrs.fa > tr_regions.fasta
+trf409.linux64 data/tr_regions.fasta 3 7 7 80 5 40 500 -h -ngs > data/grch38.tandemrepeatfinder.txt
+```
+
+Then run TRF on the reference sequence of regions:
+```bash
+trf409.linux64 data/tr_regions.fasta 3 7 7 80 5 5 500 -h -ngs > data/grch38.tandemrepeatfinder.txt
+```
+
+Because we're going to be using the variants to filter these repeat annotations, we lower the min-score from 5 to 40
+with the idea being we're more interested in sensitivity.
+

regions/intersection/variant_region_intersection.py

@@ -28,20 +28,21 @@ def main(in_bed, in_vcf, out_name):
     variants = pysam.VariantFile(in_vcf)
     with open(f"counts_{out_name}", 'w') as fout:
         for line in optional_compressed_fh(in_bed):
-            chrom, start, end = line.strip().split('\t')[:3]
+            line = line.strip()
+            chrom, start, end = line.split('\t')[:3]
             start = int(start)
             end = int(end)
             cnt = 0
             bases = 0
             for i in variants.fetch(chrom, int(start), int(end)):
                 # check only svs.. take this out for core analysis but keep in for extra analysis
-                if 'SVLEN' not in i.info or i.info["SVLEN"] < 50:
-                    continue
+                #if 'SVLEN' not in i.info or i.info["SVLEN"] < 50:
+                    #continue
                 vs, ve = truvari.entry_boundaries(i)
                 if start <= vs and ve <= end:
                     cnt += 1
                     bases += truvari.entry_size(i)
-            fout.write(f"{chrom}\t{start}\t{end}\t{cnt}\t{bases}\n")
+            fout.write(f"{line}\t{cnt}\t{bases}\n")
     data = pd.read_csv(f"counts_{out_name}", sep='\t', header=None, names=['chrom', 'start', 'end', 'num_vars', 'num_bases'])
     print("statistic\tcount\tpercent")
 

variants/README.md

@@ -106,11 +106,12 @@ First, we make the header
 bcftools merge -m none hapo_merged/*.vcf.gz -o pVCFs/GRCh38.variants.header.vcf --print-header --force-samples
 ```
 Our run had redundant sample names which were manually altered inside of the output header e.g. NA24385 is assembled by
-both eichler and li, therefor the li was renamed in the new header li:NA24385.
+both eichler and li, therefore the li was renamed in the new header li:NA24385.
 
 Next, we merge the vcfs
 ```bash
-bcftools merge -m none hapo_merged/*.vcf.gz --use-header pVCFs/GRCh38.variants.header.vcf | bcftools view -S pVCFs/sample_order.txt -o pVCFs/GRCh38.variants.vcf.gz -O z
+bcftools merge -m none hapo_merged/*.vcf.gz --use-header pVCFs/GRCh38.variants.header.vcf \
+	| bcftools view -S pVCFs/sample_order.txt -o pVCFs/GRCh38.variants.vcf.gz -O z
 ```
 The sample_order.txt is just the list of the sample names so that we can control their order in the pVCF.
 

variants/scripts/annotate_pvcf_cov.py

@@ -41,10 +41,6 @@
         if sample + '.1' not in dn_pos_coverage:
             logging.warning("missing dn %s at %s:%d.1", sample, entry.chrom, entry.start)
 
-        if entry.samples[sample]["FT"] != '.':
-            # Trust the earlier step got the coverage right
-            continue
-        
         u_cov1 = int(up_pos_coverage[sample + '.0'])
         u_cov2 = int(up_pos_coverage[sample + '.1'])
         d_cov1 = int(dn_pos_coverage[sample + '.0'])

variants/scripts/map_haplo.sh

@@ -2,20 +2,19 @@
 
 
 # Input fasta haplotype
-fasta=$1
+fasta=$(realpath $1)
 # Reference
-ref=$2
+ref=$(realpath $2)
 sample_name=$3
 # Directory to write the output folder
-out_dir=$4
+out_dir=$(realpath $4)
 
 # These would need to be a parameter / config
 params='-cx asm20 -m 10000 -z 10000,50 -r 50000,2000000 --end-bonus=100 --rmq=yes -O 5,56 -E 4,1 -B'
 threads=8
 # minimum quality score of alignments to consider
 min_qual=60
 
-
 DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
 if ! command -v minimap2 &> /dev/null
 then
@@ -25,18 +24,20 @@ fi
 
 if ! command -v paftools.js &> /dev/null
 then
-    echo "Program 'minimap2' not found in environment"
+    echo "Program 'paftools.js' not found in environment"
     exit 1
 fi
 
-
 mkdir -p ${out_dir}
 cd ${out_dir}
 anno=$DIR/annotate_cov.py
 minimap2 ${params} -t${threads} --secondary=no --cs ${ref} ${fasta} \
     | sort -k6,6 -k8,8n > aln.paf
 paftools.js stat aln.paf > aln.paf.stats.txt
-cat aln.paf | paftools.js call -f ${ref} -q ${min_qual} -L10000 - -s ${sample_name} | vcf-sort | bgzip > aln.vcf.gz
+cat aln.paf | paftools.js call -f ${ref} -q ${min_qual} -L10000 -s ${sample_name} - \
+    | vcf-sort \
+    | bcftools +fill-from-fasta /dev/stdin -- -c REF -f ${ref} \
+    | bgzip > aln.vcf.gz
 tabix aln.vcf.gz
 
 awk -v mq=${min_qual} '{if ($12 >= mq) print $6 "\t" $8 "\t" $9}' aln.paf > aln.bed

variants/scripts/mk_mergehaps.sh

@@ -1,10 +1,12 @@
+pairs=$1
 mkdir -p jobs
 mkdir -p logs
 
-cat metadata/hap_pairs.txt | while read h1 c1 h2 c2 proj samp
+DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
+cat $pairs | while read h1 c1 h2 c2 proj samp
 do
     echo "#!/bin/bash" > jobs/hapmerge_${proj}_${samp}.sh
-    echo "bash scripts/merge_haps.sh $h1 $c1 $h2 $c2 $samp hapo_merged/${proj}_${samp}.vcf.gz" >> jobs/hapmerge_${proj}_${samp}.sh
+    echo "bash $DIR/merge_haps.sh $h1 $c1 $h2 $c2 $samp hapo_merged/${proj}_${samp}.vcf.gz" >> jobs/hapmerge_${proj}_${samp}.sh
 done