Skip to content

Instantly share code, notes, and snippets.

@gedankenstuecke
Created June 17, 2015 12:56
Show Gist options
  • Save gedankenstuecke/f237566c2a154e1795b5 to your computer and use it in GitHub Desktop.
Save gedankenstuecke/f237566c2a154e1795b5 to your computer and use it in GitHub Desktop.
# this one works, plink is in cwd + conversion to BED files works flawlessly
for f in usable_files:
# gid is the genotype-ID
gid = f.split("/")[-1].split("_")[1].replace("file","")
# converts the genotyping file to plink format, using the gid as sample name
call = "./plink --23file "+ f + " F" + gid + "ID" + gid + "I 1"
call += " --out 23andme_plink/genotypeid_" + gid
print "convert gid " + gid
subprocess.call(call,shell=True)
# this is part of the same script, plink is still in cwd but it crashes, saying
# "Error: Failed to open /home/bastia.fam"
for i in list_bed:
# check that files have been processed and are working so far.
if os.path.isfile(i.replace(".bed",".fam")) and os.path.isfile(i.replace(".bed",".bim")):
listhandle.write(cwd + "/"+i.replace(".bed","")+"\n")
call_merge = "./plink --bfile " + cwd + "/" + start_bed + " --merge-list " + cwd + "/merge_list.txt --make-bed --out " + cwd + "/23andme_merged/merge-pass1"
print "merging 23andme data, first pass"
subprocess.call(["./plink","--bfile",cwd+"/"+start_bed,"--merge-list",cwd + "/merge_list.txt","--make-bed","--out",cwd+"/23andme_merged/merge-pass1"])
@audy
Copy link

audy commented Jun 17, 2015

I refactored a bit so I could read it ;)

What happens when you change line 21 to the following snippet? If you copy/paste the command does it run?

cmd = ''.join([
    './plink',
    '--bfile %s/%s' % (cwd, start_bed),
    '--merge-list %s/%s' % (cwd, '/merge_list.txt'),
    '--make-bed',
    '--out %s/%s' % (cwd, '23andme_merged/merge-pass1')
    ])

# if subprocess.call doesn't work, copy paste this to debug
print cmd

subprocess.call(cmd, shell=True)

@gedankenstuecke
Copy link
Author

Nope, the problem remains. Actually first you have to change "--bfile" etc to " --bfile" to make sure it won't concatenate into a single whitespace-free string, but even then it won't work.

@audy
Copy link

audy commented Jun 17, 2015

import subprocess
import os
from glob import glob

usable_files = glob('files/*')

# this one works, plink is in cwd + conversion to BED files works flawlessly
for f in usable_files:
    gid = os.path.basename(f).split('.')[-1]
    call = 'plink2 --23file %s F%sID%s I 1 --out 23andme_plugin/genotypeid_%s' \
            % (f, gid, gid, gid)
    subprocess.check_call(call, shell=True)

# this is part of the same script, plink is still in cwd but it crashes, saying
# "Error: Failed to open /home/bastia.fam"

cwd = os.path.dirname(os.path.abspath(__file__))
list_bed = glob('23andme_plugin/*bed')

with open('merge_list.txt', 'w') as listhandle:
    for i in list_bed:
        # check that bed and bim files exist,
        # write bed file to merge_list
        print i
        if os.path.isfile(i.replace('.bed', '.fam')) \
                and os.path.isfile(i.replace('.bed', '.bim')):
            listhandle.write(cwd + '/' + i.replace('.bed', '') + '\n')

start_bed = list_bed[0].replace('.bed', '')

cmd = ' '.join([
    'plink2',
    '--bfile %s/%s' % (cwd, start_bed),
    '--merge-list %s/%s' % (cwd, '/merge_list.txt'),
    '--make-bed',
    '--out %s/%s' % (cwd, '23andme_merged/merge-pass1')
    ])

# if subprocess.call doesn't work, copy paste this to debug
print cmd

subprocess.call(cmd, shell=True)

@audy
Copy link

audy commented Jun 17, 2015

Put the input files in files. Make sure you directories 23andme_merged and 23andme_plugin exist before running.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment