corbinq / apex

Toolkit for QTL mapping and meta-analysis.
https://corbinq.github.io/apex/
16 stars 1 forks source link

"--gene" flag in apex trans not working. #6

Closed hsun3163 closed 3 years ago

hsun3163 commented 3 years ago

Dear APEX authors,

In apex trans there is a --gene flag to include only the gene of interests.

--gene                            Restrict analysis to specified genes
                                          (gene name or comma-separated list).

However, including it doesn't seem to changes the number of genes that I analyzed, as illustrated below:

This commands supposedly only analysis one gene,

/apex trans --vcf /mnt/mfs/statgen/neuro-apex/ROSMAP-vcf/ROSMAP_chr19.vcf.gz \
--bed /mnt/mfs/statgen/neuro-apex/pipeline_testing/cache/test.chr19.mol_phe.bed.gz \
--cov /mnt/mfs/statgen/neuro-apex/pipeline_testing/cache/test.cov.gz \
--out /mnt/mfs/statgen/neuro-apex/pipeline_testing/test.19 \
--theta-file  /mnt/mfs/statgen/neuro-apex/pipeline_testing/cache/test.chr19.mol_phe.theta.gz \
--gene 'ENSG00000064666'

but instead, the function still analysis all the genes in the input bed file:

Using 48 threads.
976562 total variants on selected chromosomes.

Found 1196 samples in bcf file ... 
Found 664 samples in covariate file ... 
Found 665 samples in expression bed file ... 
Found 664 samples in common across all three files.

Processed data for 5 covariates across 664 samples.
Processed expression for 5 genes across 664 samples.
Processed genotype data for 976562 variants.

Started trans-QTL analysis ...
Scaling expression traits ... 
Calculating genotype-covariate covariance...
Calculating genotype residual variances ...Done.
Calculating expression residuals...
Scaling expression residuals ...
Processed 9766 genotype blocks out of 9765 total

Trying to use a file that contains the gene_ID to be included doesn't work either.

Could you provide an example to illustrate what is the proper way of doing the gene selection?

corbinq commented 3 years ago

Thanks for the bug report. The --gene subsetting option is not yet implemented, and should not yet be appearing in the help menu -- I will make a note to fix this.

To subset to a gene of interest, I recommend using the --bed-region argument to specify the genomic region around the gene of interest.

On Mon, Jul 19, 2021 at 2:39 PM hsun3163 @.***> wrote:

In apex trans there is a --gene flag to include only the gene of interests.

--gene Restrict analysis to specified genes (gene name or comma-separated list).

However, including it doesn't seem to changes the number of genes that I analyzed, as illustrated below:

This commands supposedly only analysis one gene,

/apex trans --vcf /mnt/mfs/statgen/neuro-apex/ROSMAP-vcf/ROSMAP_chr19.vcf.gz \ --bed /mnt/mfs/statgen/neuro-apex/pipeline_testing/cache/test.chr19.mol_phe.bed.gz \ --cov /mnt/mfs/statgen/neuro-apex/pipeline_testing/cache/test.cov.gz \ --out /mnt/mfs/statgen/neuro-apex/pipeline_testing/test.19 \ --theta-file /mnt/mfs/statgen/neuro-apex/pipeline_testing/cache/test.chr19.mol_phe.theta.gz \ --gene 'ENSG00000064666'

but instead, the function still analysis all the genes in the input bed file:

Using 48 threads. 976562 total variants on selected chromosomes.

Found 1196 samples in bcf file ... Found 664 samples in covariate file ... Found 665 samples in expression bed file ... Found 664 samples in common across all three files.

Processed data for 5 covariates across 664 samples. Processed expression for 5 genes across 664 samples. Processed genotype data for 976562 variants.

Started trans-QTL analysis ... Scaling expression traits ... Calculating genotype-covariate covariance... Calculating genotype residual variances ...Done. Calculating expression residuals... Scaling expression residuals ... Processed 9766 genotype blocks out of 9765 total

Trying to use a file that contains the gene_ID to be included doesn't work either.

Could you provide an example to illustrate what is the proper way of doing the gene selection?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/corbinq/apex/issues/6, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAWQAM6B5R3DGPANZQ3EUHDTYRWPNANCNFSM5AUH7NBQ .

hsun3163 commented 3 years ago

Thanks for the bug report. The --gene subsetting option is not yet implemented, and should not yet be appearing in the help menu -- I will make a note to fix this. To subset to a gene of interest, I recommend using the --bed-region argument to specify the genomic region around the gene of interest.

Thank you for the prompt reply, I will take a look to see what I can do.