PriceLab / TReNA-Legacy

Fit transcriptional regulatory networks using gene expression, priors, machine learning
2 stars 1 forks source link

running a function--getTfbsCountsInPromoters--in makeTRN-proximalonly-demo.R #3

Open dauss75 opened 8 years ago

dauss75 commented 8 years ago

I am interactively running the makeTRN-proximalonly-demo.R script which Seth fixed a database compatibility issue. It looks like it doesn't find any binding site for TF. I wonder whether this issue is due to the current database setup on our side.

genome.db.uri = "postgres://localhost/hg38" project.db.uri = "postgres://localhost/wholebrain" cores = 5

ptm <- proc.time()

promoter_counts = getTfbsCountsInPromoters( genome.db.uri=genome.db.uri , project.db.uri=project.db.uri , size.upstream = 10000 , size.downstream = 10000 , # define the size of the window around the TSS cores = cores ) # specify the number of cores for parallelization

proc.time() - ptm> > + + + no gene list is given, using all genes from obj@genome.db

no tf list is given. using all tfs from obj@project.db fetching promoter regions for all genes counting binding sites for each tf in each region Error in { : task 482 failed - "argument is of length zero" In addition: Warning messages: 1: closing unused connection 7 (<-localhost:11920) 2: closing unused connection 6 (<-localhost:11920) 3: closing unused connection 5 (<-localhost:11920) 4: closing unused connection 4 (<-localhost:11920) 5: closing unused connection 3 (<-localhost:11920)

user system elapsed 11.815 2.101 883.211

sessionInfo() R version 3.2.5 (2016-04-14) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 14.04.1 LTS

locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets [8] methods base

other attached packages: [1] TReNA_0.99.51 doParallel_1.0.10 iterators_1.0.8 [4] GenomicRanges_1.22.4 GenomeInfoDb_1.6.3 IRanges_2.4.8 [7] S4Vectors_0.8.11 BiocGenerics_0.16.1 RPostgreSQL_0.4-1 [10] DBI_0.4-1 vbsr_0.0.5 randomForest_4.6-12 [13] glmnet_2.0-5 foreach_1.4.3 Matrix_1.2-6

loaded via a namespace (and not attached): [1] lattice_0.20-33 codetools_0.2-14 grid_3.2.5 zlibbioc_1.16.0 [5] XVector_0.10.0 compiler_3.2.5

seth-ament commented 8 years ago

Hi Según,

Not entirely sure why it isn't working, though I have some thoughts. What do the first few rows of each of the tables look like in your version of the database?

Seth

On Mon, Jul 18, 2016 at 12:37 PM, Segun Jung notifications@github.com wrote:

I am interactively running the makeTRN-proximalonly-demo.R script which Seth fixed a database compatibility issue. It looks like it doesn't find any binding site for TF. I wonder whether this issue is due to the current database setup on our side.

genome.db.uri = "postgres://localhost/hg38" project.db.uri = "postgres://localhost/wholebrain" cores = 5

ptm <- proc.time()

promoter_counts = getTfbsCountsInPromoters( genome.db.uri=genome.db.uri , project.db.uri=project.db.uri , size.upstream = 10000 , size.downstream = 10000 , # define the size of the window around the TSS cores = cores ) # specify the number of cores for parallelization

proc.time() - ptm> > + + + no gene list is given, using all genes from obj@genome.db

no tf list is given. using all tfs from obj@project.db fetching promoter regions for all genes counting binding sites for each tf in each region Error in { : task 482 failed - "argument is of length zero" In addition: Warning messages: 1: closing unused connection 7 (<-localhost:11920) 2: closing unused connection 6 (<-localhost:11920) 3: closing unused connection 5 (<-localhost:11920) 4: closing unused connection 4 (<-localhost:11920) 5: closing unused connection 3 (<-localhost:11920)

user system elapsed 11.815 2.101 883.211

sessionInfo() R version 3.2.5 (2016-04-14) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 14.04.1 LTS

locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets [8] methods base

other attached packages: [1] TReNA_0.99.51 doParallel_1.0.10 iterators_1.0.8 [4] GenomicRanges_1.22.4 GenomeInfoDb_1.6.3 IRanges_2.4.8 [7] S4Vectors_0.8.11 BiocGenerics_0.16.1 RPostgreSQL_0.4-1 [10] DBI_0.4-1 vbsr_0.0.5 randomForest_4.6-12 [13] glmnet_2.0-5 foreach_1.4.3 Matrix_1.2-6

loaded via a namespace (and not attached): [1] lattice_0.20-33 codetools_0.2-14 grid_3.2.5 zlibbioc_1.16.0 [5] XVector_0.10.0 compiler_3.2.5

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/PriceLab/TReNA/issues/3, or mute the thread https://github.com/notifications/unsubscribe-auth/AG7ejpc1zcBlhVPoqgMnO0n3jlXN4wXZks5qW6uzgaJpZM4JO5vv .

Seth A. Ament, Ph.D. | Postdoctoral Fellow Institute for Systems Biology | 401 Terry Ave N, Seattle, WA 98109 p. 206.732.1346 www.systemsbiology.org

Revolutionizing science. Enhancing life.

This message (including any attachments) contains confidential information intended for a specific individual and purpose, and is protected by law. If you are not the intended recipient, you should delete this message and are hereby notified that any disclosure, copying, or distribution of this message, or the taking of any action based on it, is strictly prohibited.

dauss75 commented 8 years ago

For wholebrain:

  1. footprints table

chr | mfpstart | mfpend | strand | motifname | motiflength | footprintlength | score | pval | qval | sequence -------+-----------+-----------+--------+-------------------------+-------------+-----------------+-------+----------+----------+------------------------------- chr4 | 96062473 | 96062491 | + | MA0681.1 | 11 | 18 | 54.3 | 3.67e-06 | 0.413 | TAACCTAATTA chr4 | 96062473 | 96062491 | + | MA0713.1 | 11 | 18 | 55.7 | 2.67e-06 | 0.402 | TAACCTAATTA chr4 | 96062473 | 96062493 | + | MA0892.1 | 10 | 18 | 51.5 | 7.08e-06 | 0.909 | CCTAATTAGA chr4 | 96062473 | 96062496 | - | MA0713.1 | 11 | 18 | 51.5 | 7.03e-06 | 0.409 | AAATCTAATTA

  1. motifsgenes

       motif             | tf_name  |     tf_ensg

    -------------------------------+----------+----------------- MA0002.2 | RUNX1 | ENSG00000159216 MA0003.3 | TFAP2A | ENSG00000137203 MA0004.1 | ARNT | ENSG00000143437 MA0006.1 | AHR | ENSG00000106546 MA0006.1 | ARNT | ENSG00000143437 MA0007.3 | AR | ENSG00000169083

For hg38:

  1. gtf

    chr | start | endpos | score | strand | frame | moleculetype | gene_id | gene_version | gene_name | gene_source | gene_biotype | havana_gene | havana_gene_version | transcript_id | transcript_version | transcript_name | transcript_source | transcript_biotype | havana_transcript | havana_transcript_version | tag | transcript_support_level | exon_number | exon_id | exon_version | ccds_id | protein_id | protein_version | annotation -------+-----------+-----------+-------+--------+-------+-----------------+-----------------+--------------+------------------------+----------------+------------------------------------+--------------------+---------------------+-----------------+--------------------+----------------------------+-------------------+------------------------------------+--------------------+---------------------------+---------------+--------------------------+-------------+-----------------+--------------+-----------+-----------------+-----------------+---------------- chr1 | 11869 | 14409 | . | + | . | gene | ENSG00000223972 | 5 | DDX11L1 | havana | transcribed_unprocessed_pseudogene | OTTHUMG00000000961 | 2 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | havana chr1 | 11869 | 14409 | . | + | . | transcript | ENSG00000223972 | 5 | DDX11L1 | havana | transcribed_unprocessed_pseudogene | OTTHUMG00000000961 | 2 | ENST00000456328 | 2 | DDX11L1-002 | havana | processed_transcript | OTTHUMT00000362751 | 1 | basic | 1 | NA | NA | NA | NA | NA | NA | havana

paul-shannon commented 8 years ago

Hi Segun,

Sorry you are having trouble. Let's figure it out...

Seth has written some nice unit tests to accompany his new runTReNA functions. Running those tests is a great place to start, specifically

test_getTfbsCountsInPromoters()

should run quickly and correctly.

The next step we will try is to make direct calls to your database, following the style and syntax Seth uses in his function.

The step after that is pending, waiting on the completion of the piq footprints. Once we have three sorts of footprints (wellington, HINT, piq) we will update the database schema, consolidate the 3 sources of footprints into one large table (if possible).

Let us know if the test function listed above runs successfully.

dauss75 commented 8 years ago

Hi Paul,

The function--test_getTfbsCountsInPromoters()-- worked fine.

seth-ament commented 8 years ago

Hi Según,

I'm really sorry this has not worked correctly. It's a top priority for me to get this working for you!

Would it be possible for you to give me a login so that I can try running TReNA functions in your environment and using your database instances?

Thanks!

Seth

On Wed, Jul 20, 2016 at 7:55 AM, Segun Jung notifications@github.com wrote:

Hi Paul,

The function--test_getTfbsCountsInPromoters()-- worked fine.

  • Segun

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/PriceLab/TReNA/issues/3#issuecomment-233974642, or mute the thread https://github.com/notifications/unsubscribe-auth/AG7ejodYZeMQ3Y6DYrUHOZCssTtQgkKVks5qXjbVgaJpZM4JO5vv .

Seth A. Ament, Ph.D. | Postdoctoral Fellow Institute for Systems Biology | 401 Terry Ave N, Seattle, WA 98109 p. 206.732.1346 www.systemsbiology.org

Revolutionizing science. Enhancing life.

This message (including any attachments) contains confidential information intended for a specific individual and purpose, and is protected by law. If you are not the intended recipient, you should delete this message and are hereby notified that any disclosure, copying, or distribution of this message, or the taking of any action based on it, is strictly prohibited.

dauss75 commented 8 years ago

I am trying to have better understanding on the overall workflow. Could you check/explain my notes below that are based on the "systems_genetics_pipeline_2016-3-11" ppt file?

1) we have tissue-specific (lymphoblast for now) candidate TF binding sites and TReNA (gene) model that will be used for "physical" TF-target gene network. Is this where we are at now?

2) I am looking into the runTReNA-demo.R script to have an idea of GUI. From the script, what I understand in terms of inputs/outputs is: a) the main function (runTReNA) requires gene, mtx and candidate.tfs from which gene is a list of genes from user; mtx is gene expression data (is this the gene model from the makeTrnFromPromoterCountsAndExpression function?); candidate.tfs is merely a list of TFs.

3) in terms of GUI design, we will need a) input: gene list from user b) output: plain text file? c) options:

Could you comment what I am missing? I am pretty sure there are plenty.

-Segun

seth-ament commented 8 years ago

Hi Según,

1) we have tissue-specific (lymphoblast for now) candidate TF binding sites and TReNA (gene) model that will be used for "physical" TF-target gene network. Is this where we are at now?

*Actually, the TReNA model we've created already incorporates information both about TF binding sites and gene co-expression.

getTfbsCountsInPromoters() takes the results of the genomic footprinting (together with gene models) and counts the number of predicted binding sites for each TF within 10kb of each transcription start site. This is the "physical" network.

makeTrnFromPromoterCountsAndExpression() performs the integration with tissue-specific expression data. In the example script, we are using gene expression microarrays from GSE37772, which are loaded as part of the TReNA package. Moving forward, we'll probably use data from GTEx, and Cory and I have used tens of other datasets in the past. For each gene, we consider a set of TFs that have strong regulatory potential at that TSS, based on the previous steps. Specifically, we rank potential target genes for each TF by the TFBS counts from step 1, and we choose the upper quartile as candidate regulators (other thresholds are possible). We then use the expression data to select a subset of these TFs as active regulators by building a model to predict the expression of the target gene based on the expression of the TFs. By default, elastic net regression is used, but the RandomForest and bayesSpike methods are also implemented. This is the final, genome-scale TRN model.

2) I am looking into the runTReNA-demo.R script to have an idea of GUI. From the script, what I understand in terms of inputs/outputs is: a) the main function (runTReNA) requires gene, mtx and candidate.tfs from which gene is a list of genes from user; mtx is gene expression data (is this the gene model from the makeTrnFromPromoterCountsAndExpression function?); candidate.tfs is merely a list of TFs

*I really like the idea of creating a GUI. Honestly, though, I wouldn't use the runTReNA-demo.R script. I think we should focus on building genome-scale models, as in makeTRN-proximalonly-demo.R. Consider as options all the variables in getTfbsCountsInPromoters() and makeTrnFromPromoterCountsAndExpression().

Hopefully this clarifies things a bit. Let's keep talking about it during our conference call tomorrow.

Seth

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/PriceLab/TReNA/issues/3#issuecomment-235693175, or mute the thread https://github.com/notifications/unsubscribe-auth/AG7ejkqGpHORmfoam6B1tZyvbSMGL-Omks5qZ7D8gaJpZM4JO5vv .

Seth A. Ament, Ph.D. | Postdoctoral Fellow Institute for Systems Biology | 401 Terry Ave N, Seattle, WA 98109 p. 206.732.1346 www.systemsbiology.org

Revolutionizing science. Enhancing life.

This message (including any attachments) contains confidential information intended for a specific individual and purpose, and is protected by law. If you are not the intended recipient, you should delete this message and are hereby notified that any disclosure, copying, or distribution of this message, or the taking of any action based on it, is strictly prohibited.

dauss75 commented 8 years ago

Hi Seth,

I am trying to incorporate the output from Cory's workflow:

money shot

bedtools intersect -wo -a /scratch/resources/GRCh38.84.1mb.bed -b /scratch/merge/all_redux2.bed > /scratch/merge/tfdb_lympho_1mb_GRCh38.84.bed

to your function: getTfbsCountsInPromoters

to create a wrapper.

By looking at the function, it doesn't seem to work directly with the output which is in bed format.

getTfbsCountsInPromoters function (genome.db.uri, project.db.uri, genelist = NULL, tflist = NULL, biotype = "protein_coding", moleculetype = "gene", use_gene_ids = T, size.upstream = 10000, size.downstream = 10000, cores = 1, verbose = 1)

I guess Cory generated the final output and passed on to you so I wonder whether you have run the getTfbsCountsInPromoters function with the input. It would be very useful if you can let me know what is the exact inputs.

Thanks,

Segun

seth-ament commented 8 years ago

Hi Según,

getTfbsCountsInPromoters() is designed to pull data from a footprints file that is incorporated into a Postgres database. This avoids having to load the entire footprint file into memory.

I believe Paul has previously sent you this URL, providing instructions for filling the database with tables, including the footprints bed file: https://github.com/PriceLab/snpFoot/tree/master/inst/misc/postgres I think these instructions should work to incorporate the new footprints file that you have generated so that you could run the entire workflow from end to end.

However, we are currently updating the schema for the footprints table so that the same fields apply to Wellington, HINT, PIQ, or an ensemble of scores from those (and other) methods. So if you wait a couple days we might be ready to hand over a new stable version of that workflow.

Seth

On Tue, Aug 2, 2016 at 1:25 PM, Segun Jung notifications@github.com wrote:

Hi Seth,

I am trying to incorporate the output from Cory's workflow: money shot

bedtools intersect -wo -a /scratch/resources/GRCh38.84.1mb.bed -b /scratch/merge/all_redux2.bed > /scratch/merge/tfdb_lympho_1mb_GRCh38.84.bed

to your function: getTfbsCountsInPromoters

to create a wrapper.

By looking at the function, it doesn't seem to work directly with the output which is in bed format.

getTfbsCountsInPromoters function (genome.db.uri, project.db.uri, genelist = NULL, tflist = NULL, biotype = "protein_coding", moleculetype = "gene", use_gene_ids = T, size.upstream = 10000, size.downstream = 10000, cores = 1, verbose = 1)

I guess Cory generated the final output and passed on to you so I wonder whether you have run the getTfbsCountsInPromoters function with the input. It would be very useful if you can let me know what is the exact inputs.

Thanks,

Segun

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/PriceLab/TReNA/issues/3#issuecomment-237032219, or mute the thread https://github.com/notifications/unsubscribe-auth/AG7ejv6UVKMJk5WaSlyV388jaDl1m1oPks5qb6esgaJpZM4JO5vv .

Seth A. Ament, Ph.D. | Postdoctoral Fellow Institute for Systems Biology | 401 Terry Ave N, Seattle, WA 98109 p. 206.732.1346 www.systemsbiology.org

Revolutionizing science. Enhancing life.

This message (including any attachments) contains confidential information intended for a specific individual and purpose, and is protected by law. If you are not the intended recipient, you should delete this message and are hereby notified that any disclosure, copying, or distribution of this message, or the taking of any action based on it, is strictly prohibited.

dauss75 commented 8 years ago

Thanks Seth for reminding me this. I haven't linked this to Paul's instruction.

Segun