import os
sdir = '/hapmap/'
basebedname = 'hapmap3'
ddir = '/erlotinib/'
temponame = 'tempo'
phenoname = 'erlotinib'
os.system('plink --bfile ' + sdir + basebedname + ' --recode12 --out ' + ddir + temponame )
## join tempo.ped with phenodata.txt
## and generate pheno.ped file
shortfilename = ddir + "phenodata.txt" ## FID,IID,dad,mom,sex,pheno1,...
longfilename = ddir + temponame + ".ped" ## FID,IID,dad,mom,sex,pheno,snp1,...
outfilename = ddir + phenoname + ".ped"
index = {}
shortfile= open( shortfilename, "r" )
for row in shortfile:
part = row.split(None,2)
iid = part[1] ## IID is the second column in phenofile
index[iid] = row.rstrip()
## rstrip gets rid of \n and other white spaces
shortfile.close()
outfile = open( outfilename, "w")
longfile= open( longfilename, "r" )
for row in longfile:
part = row.split(None,6)
iid = part[1]
restofline = part[6]
if iid in index:
pheno = index[iid]
outfile.write( pheno + ' ' + restofline)
longfile.close()
outfile.close()
## generate pheno.dat file
infilename = ddir + temponame + '.map'
outfilename = ddir + phenoname + '.dat'
infile = open(infilename,'r')
outfile = open(outfilename,'w')
## write covariates and phenotype names
outfile.write('C igrowth\n')
outfile.write('T logAUC\n')
## write one row per snp
for linea in infile:
part = linea.split(None,2)
newlinea = 'M ' + part[1] + '\n'
outfile.write(newlinea)
infile.close()
outfile.close()
Friday, April 20, 2012
python script to convert plink bed/bim to qtdt format
Friday, April 13, 2012
create new theme for ggplot
from: http://sape.inf.usi.ch/quick-reference/ggplot2/themes
theme_complete_bw <- function(base_size = 12) {
structure(list(
axis.line = theme_blank(),
axis.text.x = theme_text(size = base_size * 0.8 , lineheight = 0.9, colour = "black", vjust = 1),
axis.text.y = theme_text(size = base_size * 0.8, lineheight = 0.9, colour = "black", hjust = 1),
axis.ticks = theme_segment(colour = "black"),
axis.title.x = theme_text(size = base_size, vjust = 0.5),
axis.title.y = theme_text(size = base_size, angle = 90, vjust = 0.5),
axis.ticks.length = unit(0.15, "cm"),
axis.ticks.margin = unit(0.1, "cm"),
legend.background = theme_rect(colour=NA),
legend.key = theme_rect(fill = NA, colour = "black", size = 0.25),
legend.key.size = unit(1.2, "lines"),
legend.text = theme_text(size = base_size * 0.8),
legend.title = theme_text(size = base_size * 0.8, face = "bold", hjust = 0),
legend.position = "right",
panel.background = theme_rect(fill = NA, colour = "black", size = 0.25),
panel.border = theme_blank(),
panel.grid.major = theme_line(colour = "black", size = 0.05),
panel.grid.minor = theme_line(colour = "black", size = 0.05),
panel.margin = unit(0.25, "lines"),
strip.background = theme_rect(fill = NA, colour = NA),
strip.text.x = theme_text(colour = "black", size = base_size * 0.8),
strip.text.y = theme_text(colour = "black", size = base_size * 0.8, angle = -90),
plot.background = theme_rect(colour = NA, fill = "white"),
plot.title = theme_text(size = base_size * 1.2),
plot.margin = unit(c(1, 1, 0.5, 0.5), "lines")
), class = "options")
}
Thursday, April 12, 2012
# cdf: cd's to frontmost window of Finder (by Yarko)
copy the following in ~/.bash_profile
# cdf: cd's to frontmost window of Finder
cdf ()
{
currFolderPath=$( /usr/bin/osascript <<" EOT"
tell application "Finder"
try
set currFolder to (folder of the front window as alias)
on error
set currFolder to (path to desktop folder as alias)
end try
POSIX path of currFolder
end tell
EOT
)
# ignore the cd if it's my desktop
if [ "$currFolderPath" != "$HOME/Desktop/" ]
then
echo "cd to \"$currFolderPath\""
cd "$currFolderPath"
else
echo "no finder path to cd to..."
fi
}
Subscribe to:
Comments (Atom)