hakyimlab/predictdbpipeline

Deprecation warning

This code has been deprecated and exists solely for reference and documentation.

The model prediction pipeline (with nested cross validation) is available in github here, and a brief overview is available here (PrediXcan google group).

GTEx Elastic-Net Model Pipeline

This directory can be used for creating models for PredictDB with GTEx data.

These models can be used to predict gene expression in 44 different tissues, based on individual genotype. To train the models, we used version 6p of the GTEx data and applied 10-fold cross validation elastic-net within SNPs that are located within 1mb of the gene as features.

The pipeline takes 4 different sets of inputs (all in the data/inputs/ directory, possibly in the form of symbolic links), does some preprocessing, stores the results in the data/intermediate/ directory. Once this is all done, we can apply our machine learning technique, output the results in the data/output/ directory. There's also a script to turn this data into SQLite databases, to be used by PrediXcan.

With some modifications, this pipeline can be extended to analyze other studies which have genotype and expression data.

See our wiki for a tutorial with example files and a detailed description of these scripts.

predictdbpipeline's Issues

Error in smooth.spline during Filtering on significance

Hello,

I got some error during runing post-process.py. This following is the error:
$ ./post_process.py
Generating ../../data/output/dbs/Test_Test_alpha0.5_window1e6.db...
Creating weights table...
Creating extra table...
Creating construction table...
Creating meta_data table...
Commiting...
Done.
Filtering Test on significance.
Error in smooth.spline(lambda, pi0, df = smooth.df) :
missing or infinite values in inputs are not allowed
Calls: package_db -> qvalue -> pi0est -> smooth.spline
In addition: Warning messages:
1: In min(p) : no non-missing arguments to min; returning Inf
2: In max(p) : no non-missing arguments to max; returning -Inf
3: In min(p) : no non-missing arguments to min; returning Inf
4: In max(p) : no non-missing arguments to max; returning -Inf
Execution halted

It look like the weight table in unfilter SQLiteDB have blank p.value and other 4 last values as well.

Thanks,
Phuwanat

error when generating database with all of the prediction models

Hi dear hakyimlab group,

I was trying to run script make_sqlite_db.py for making database with all of the prediction models produced by glmnet. And, I got this error:

[rchen7@hii example2]$ python ../../scripts/make_sqlite_db.py --output ../../data/output/dbs/last_tw_elasticnet.db --betas ../../data/output/allBetas/allbetas --results ../../data/output/allResults/allresults --construction ../../data/output/allLogs/alllogs --meta ../../data/output/allMetaData/last.allMetaData.txt

Generating ../../data/output/dbs/last_tw_elasticnet.db...
Creating weights table...
Traceback (most recent call last):
File "../../scripts/make_sqlite_db.py", line 187, in
args.meta, args.output, args.intron)
File "../../scripts/make_sqlite_db.py", line 143, in make_sqlite_db
db.add_weights(betas)
File "../../scripts/make_sqlite_db.py", line 34, in add_weights
self.insert_weights_row(row)
File "../../scripts/make_sqlite_db.py", line 39, in insert_weights_row
(row['rsid'], row['gene'], row['beta'], row['ref'],row['alt']))
KeyError: 'rsid'

Could anyone advice on this? Really appreciate!!

Ruoxi

convert genotype file

Hi,
I have a genotype file in VCF 4.1 format (get from GATK), can I convert this file to the INPUT genotype file simply recoding as below,
0/0 ---->0;
0/1 ---->1;
1/1 ---->2;
Thank you!

code to generate SNP annotation file

Is there any code to generate the SNP annotation file? Thanks!

Best,
Ya

train_models issue

Hi!
I am going through the tutorial right now but when I'm running the train _models step, it just create empty geuvadis files .risoe-r12-cn534 and I don't think this is the outcome expected since I can't pursue with the rest of the pipeline.
Do you know where that problem could came from? I am expecting maybe an error in the qsub command but since it's still producing the files in the end I'm not sure...

Thank you!

command qsub not found

Hello,
The computing cluster I am working with uses Slurm instead of SGE and requires the $SBATCH command to submit jobs to the queue. Because of this I am running into an error when attempting to run the train_models.py file. Is there a simple way to alter this file so I can run it on a slurm based system?

empty value in genotype file

Hi dear hakyimlab, @hakyim
In my genotype file, there are some value are empty, so i use "NA" to replace them, like the following . Does this matter?
4_165621955_A_C_b37 0 NA 0 0 0 0 NA 1 1 0 0 0 0 0 0 NA 0 0 0 0 0 0 NA 0 0 0 0 NA 0 0 0 0 0 0 0 NA 0 0 0 0 1 1 0 0 0 0 NA NA 1 0 1 0 0 0 NA 0 0 NA NA 0 1 0 0 0 0 0 2 0 1 0 0 0 0 0 0 NA 1 0 NA 0 0 NA NA 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 NA NA 0 0 0 0 0 NA 0 1 0 1 0 0 0 0 2 0 0 0 0 0 0 0 NA 1 0 0 0 0 0 1 0 0 0