search

LieberInstitute / zandiHyde_bipolar_rnaseq

RNA-seq analysis for psychENCODE R01
6 stars 1 forks source link

PGC TWAS results into a single Rdata file #2

Closed lcolladotor closed 3 years ago

lcolladotor commented 3 years ago

Hi Arta,

Can you adapt https://github.com/LieberInstitute/twas/blob/master/bsp2/read_twas.R#L19-L228 to combine into a single Rdata file all the results from PGC's BPD GWAS results from our TWAS pipeline applied to sACC and Amygdala? Thanks!

Best, Leo

PS I'm trying to make some issues to list some of the steps remaining.

aseyedia commented 3 years ago

Made some progress on this; Just for clarification, which rse object should I use for the following?

https://github.com/LieberInstitute/twas/blob/c9a4734b685c36c20a92682ce1b2b7e0158b7ba0/bsp2/read_twas.R#L149-L153

We have multiple rse_genes: one source rse_gene for each brain region, amygdala and sacc, and the one input into filter_snps.R:

https://github.com/LieberInstitute/zandiHyde_bipolar_rnaseq/blob/57427e1c0f15dadf891176badc3d0b89fe7fddc9/dev_twas/filter_data/filter_snps.R#L56-L66

and another that is saved after being modified:

https://github.com/LieberInstitute/zandiHyde_bipolar_rnaseq/blob/57427e1c0f15dadf891176badc3d0b89fe7fddc9/dev_twas/filter_data/filter_snps.R#L164-L175

I'm guess I'm just making sure which rse_gene object I should use and asking if I should save everything as a single Rdata file for both sacc and amygdala, or if there should be two separate ones, one for each brain subregion.

lcolladotor commented 3 years ago

Hi Arta,

You can use the sacc ones if you want, or the amygdala ones. It doesn't matter. The RSE files in this script are really only used for adding the gencode gene ID and the gene symbol as in https://github.com/LieberInstitute/twas/blob/master/bsp2/read_twas.R#L191-L220.

Having said that, it's maybe easier to load the subsetted RSE files that were used in the build_bims script. The only one that might be a bit harder is the exon-exon junction RSE in case the sACC and Amygdala have different exon-exon junctions. I forget if Emily had already matched them or not.

Best, Leo

aseyedia commented 3 years ago

I didn't fundamentally change anything about the script, but I did add flag options to the script with an automatic assignment of sacc to opt$region just in case it mattered. The Rdata file is stored in /dcl01/lieber/ajaffe/lab/zandiHyde_bipolar_rnaseq/dev_twas/rda/twas_exp.Rdata.