Shapeit arises an error for provided pipline

[sinamajidian]

Hi
My question may be out of scope of your program.
In the pipeline you mentioned
shapeit --input-vcf VCFfile -M genetic_map.txt --output-graph VCFfile.graph
I faced with this error
ERROR: You cannot phase less than 10 samples without using a reference panel!
Should I add --input-ref ?

[vibansal]

Yes, you have to provide the reference panel files for running shapeit. I have updated the instructions.

[sinamajidian]

Thanks. Now, I'm trying run step 3 of the pipeline python samplehaps.py VCFfile

File "../IntegratedPhasing-master/samplehaps.py", line 29
samplehap = SHAPEIT + ' -convert' + ' --seed ' + `r` + ' --input-graph ' + VCF + '.graph' + ' --output-sample ' + outfile + ' -L ' + logfile;
                                                     ^
SyntaxError: invalid syntax

it seems that it should be 'r' rather than `r ` . After fixing that, when I run python samplehaps.py VCFfile The output is

program to sample haplotypes from SHAPEIT2 HMM 
set path to shapeit before running this program, the shapeit graph file also needs to exist 
usage: python samplehaps.py VCFfile number_haps(integer)

I think it needs three arguments as it mentioned in line 46

if len(sys.argv) < 3:

    sys.exit();
sample_haplotypes(sys.argv[1],int(sys.argv[2])); 

I think its nsamples as it called in line 7 def sample_haplotypes(VCF,nsamples):
Is it related to 1000g dataset? I used 10 for now.

$ python ../IntegratedPhasing-master/samplehaps.py VCFfile 10
Traceback (most recent call last):
  File "../IntegratedPhasing-master/samplehaps.py", line 53, in <module>
    sample_haplotypes(sys.argv[1],int(sys.argv[2])); 
  File "../IntegratedPhasing-master/samplehaps.py", line 13, in sample_haplotypes
    VCF_file = open(VCF);
IOError: [Errno 2] No such file or directory: 'VCFfile'

I think I should changed VCF_file = open(VCF); to VCF_file = open(VCF+'.vcf');
Then i faced wit NameError: name 'xrange' is not defined. I think I should use python2

Converter mode : the argument ('r') for option '--seed' is invalid
Traceback (most recent call last):
  File "../IntegratedPhasing-master/samplehaps.py", line 53, in <module>
    sample_haplotypes(sys.argv[1],int(sys.argv[2])); 
  File "../IntegratedPhasing-master/samplehaps.py", line 32, in sample_haplotypes
    File = open(outfile + '.haps','r');
IOError: [Errno 2] No such file or directory: 'VCFfile.phased.haps'

It needs 'VCFfile.phased.haps' . Would you please tell me how I should provide this file?

Ps. I think +`r`+ problem was for python2. I'm familiar with Python3 version in which +str(r)+ is used. Now, it's working. But output is empty. Maybe my vcf file has problem.

Now, I'm using vcf provided in DATA folder. After running third step python samplehaps.py VCFfile 1000 > VCFfile.haps . When I check VCFfile.haps

(base) majid001@waterman:/mnt/LTR_userdata/majid001/a/b$ head -n 23 VCFfile.haps

Segmented HAPlotype Estimation & Imputation Tool
  * Authors : Olivier Delaneau, Jared O'Connell, Jean-François Zagury, Jonathan Marchini
  * Contact : send an email to the OXSTATGEN mail list https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=OXSTATGEN
  * Webpage : https://mathgen.stats.ox.ac.uk/shapeit
  * Version : v2.r904
  * Date    : 08/08/2019 16:59:16
  * LOGfile : [VCFfile.shapeit.log]

MODE -convert : EXTRACTING HAPLOTYPES FROM GRAPH FILES

Parameters :
  * Seed : 0
  * Parallelisation: 1 threads
  * Ref allele is NOT aligned on the reference genome

Reading genotype graphs in [VCFfile.graph]

ERROR: Binary graph file format unsupported!

Oh! VCFfile.graph is empty. Second step (shapeit) was not completed successfully.

ERROR: Reference and Main panels are not well aligned:
  * #Missing sites in reference panel = 989
  * #Misaligned sites between panels = 100
  * #Multiple alignments between panels = 0

I think I should use --exclude-snp *.snp.strand.exclude.

[sinamajidian]

Output of third step is two files VCFfile.haps

Segmented HAPlotype Estimation & Imputation Tool
  * Authors : Olivier Delaneau, Jared O'Connell, Jean-François Zagury, Jonathan Marchini
  * Contact : send an email to the OXSTATGEN mail list https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=OXSTATGEN
  * Webpage : https://mathgen.stats.ox.ac.uk/shapeit
  * Version : v2.r904
  * Date    : 08/08/2019 17:22:23
  * LOGfile : [VCFfile.shapeit.log]
MODE -convert : EXTRACTING HAPLOTYPES FROM GRAPH FILES
Parameters :
  * Seed : 0
  * Parallelisation: 1 threads
  * Ref allele is NOT aligned on the reference genome

and VCFfile.hapsamples

2 3 4 8 9 10 11 12 13 15 18 
0100000001100000000010000
0100000001100000000010000

When I run the fourth step, I faced this

$ python2 ../IntegratedPhasing-master/encodereads.py VCFfile.haps > sample.pseudo_reads

read 2999 haplotypes from file
Traceback (most recent call last):
  File "../IntegratedPhasing-master/encodereads.py", line 74, in <module>
    else: summarize_haps(haplist,varindex);
  File "../IntegratedPhasing-master/encodereads.py", line 24, in summarize_haps
    counts = pairwiseLD(haplist,nsamples,i,j);
  File "../IntegratedPhasing-master/encodereads.py", line 7, in pairwiseLD
    if haplist[r][i] == '0' and haplist[r][j] == '0': counts[0] +=1; 
IndexError: string index out of range

It seems that it the argument should be VCFfile.hapsamples rather than VCFfile.haps.

Yes, It works!

head -n 10 sample.pseudo_reads 
2 FRAGMENT_1 0 -1 -1 2 0 3 1 ;;
2 FRAGMENT_2 0 -1 -1 2 0 4 0 ;;
2 FRAGMENT_3 0 -1 -1 2 0 8 0 ;;
2 FRAGMENT_4 0 -1 -1 2 0 9 0 ;;

To summarize,
1- change in samplehaps.py VCF_file = open(VCF+'.vcf')
2- Run step 2 twice. For the second time, --exclude-snp *.snp.strand.exclude should be added.
3- step 4 should be changed to encodereads.py VCFfile.hapsamples > sample.pseudo_reads
and maybe in step 3 python samplehaps.py VCFfile 1000 is enough. > VCFfile.haps is not needed or maybe changed to > VCFfile.haps.log

[vibansal]

Thanks for debugging the pipeline. I have updated the instructions.

[sinamajidian]

You're welcome.
When > VCFfile.haps is removed in step 3, it becomes slow. It seems that the code needs less time to run rather than print on screen some lines for each sample. so maybe it's better to use > VCFfile.log

To make the code samplehaps.py campatible with python3 +`r` + in line 27 should be replaced with +str(r)+ .

Maybe helpful for one another who wants to use the pipeline:
I faced another error using shapeit

Checking missingness and MAF...
ERROR: 20 SNPs with high rates of missing data (>10%).  These sites should be removed.

I found out that the problem was with the non-genotyped variants in the VCF file. so I removed them using grep -v "\.:\.:" all.vcf > genotyped.vcf. Now it's working!