Same as 8b53f1e but done differently

- support multi-sample VCFs but at the plotting stage, the roh detection stage was fine
samtools · Sep 8, 2017 · d0f5845 · d0f5845
1 parent 8b53f1e
commit d0f5845
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 35 deletions.
diff --git a/misc/plot-roh.py b/misc/plot-roh.py
@@ -263,22 +263,24 @@ def parse_samples(fname,highlight):
             pos  = int(row[2])
             reg  = region_overlap(regs,chr,pos,pos)
             if reg==None: continue
-            smpl = row[3]
-            if samples!=None and smpl not in samples: continue
-            gt   = row[4]
-            x = gt.split('/')
-            dsg = 2
-            if x[0]!=x[1]: dsg = 1
-            elif x[0]=='0': continue        # skip HomRef 0/0 genotypes
-            if chr not in dat_gt: 
-                dat_gt[chr] = {}
-                chrs.append(chr)
-            if smpl not in dat_gt[chr]: 
-                dat_gt[chr][smpl] = []
-            if smpl not in smpl2y:
-                y = len(smpl2y)
-                smpl2y[smpl] = y
-            dat_gt[chr][smpl].append([pos,dsg])
+            for i in range(3,len(row),2):
+                smpl = row[i]
+                if samples!=None and smpl not in samples: continue
+                gt   = row[i+1]
+                x = gt.split('/')
+                if x[0]=='.': continue          # missing genotype ./.
+                dsg = 2
+                if x[0]!=x[1]: dsg = 1
+                elif x[0]=='0': continue        # skip HomRef 0/0 genotypes
+                if chr not in dat_gt: 
+                    dat_gt[chr] = {}
+                    chrs.append(chr)
+                if smpl not in dat_gt[chr]: 
+                    dat_gt[chr][smpl] = []
+                if smpl not in smpl2y:
+                    y = len(smpl2y)
+                    smpl2y[smpl] = y
+                dat_gt[chr][smpl].append([pos,dsg])
         elif row[0]=='RG':
             smpl  = row[1]
             if samples!=None and smpl not in samples: continue
@@ -326,8 +328,8 @@ def parse_samples(fname,highlight):
     off += max_pos + off_sep
     off_list.append(off)
 
-height = len(fnames)
-if len(fnames)>5: heigth = 5
+height = len(smpl2y)
+if len(smpl2y)>5: heigth = 5
 wh = 20,height
 
 def bignum(num):

diff --git a/misc/run-roh.pl b/misc/run-roh.pl
@@ -172,23 +172,13 @@ sub run_roh
     while (my $file = readdir($dh))
     {
         if ( !($file=~/\.vcf$/i) && !($file=~/\.vcf\.gz$/i) && !($file=~/\.bcf$/i) ) { next; }
-        my $bname = $`;
-        my @samples = cmd("bcftools query -l '$$opts{indir}/$file'");
-        for my $sample (@samples)
-        {
-            $sample =~ s/\s*$//;
-            my $smpl = $sample;
-            $smpl =~ s{/+}{_}g;
-            $smpl =~ s{\s+}{_}g;
-            my $outfile = "$$opts{outdir}/$bname.$smpl.bcf";
-            push @files,$outfile;
-            if ( -e $outfile ) { next; }
-            cmd(
-                "bcftools view -s '$sample' '$$opts{indir}/$file' -Ou | " .
-                "bcftools annotate --rename-chrs $chr_fname -Ou | " .
-                "bcftools annotate -c CHROM,POS,REF,ALT,AF1KG -h $$opts{af_annots}.hdr -a $$opts{af_annots} -Ob -o $outfile.part && " .
-                "mv $outfile.part $outfile",%$opts);
-        }
+        my $outfile = "$$opts{outdir}/$`.bcf";
+        push @files,$outfile;
+        if ( -e $outfile ) { next; }
+        cmd(
+            "bcftools annotate --rename-chrs $chr_fname '$$opts{indir}/$file' -Ou | " .
+            "bcftools annotate -c CHROM,POS,REF,ALT,AF1KG -h $$opts{af_annots}.hdr -a $$opts{af_annots} -Ob -o $outfile.part && " .
+            "mv $outfile.part $outfile",%$opts);
     }
     closedir($dh) or error("close failed: $$opts{indir}");