cannin/biostars_qa huggingface dataset:
Question: I have installed GenomicRanges, the dependency before-hand and then have tried to install GenomicFeatures. Both appear to say that they have successfully been installed however RStudio does not recognise any of them.

Below are the commands and outputs I have used:

> ## try http:// if https:// URLs are not supported
> source("https://bioconductor.org/biocLite.R")
Bioconductor version 3.6 (BiocInstaller 1.28.0), ?biocLite for help
> biocLite("GenomicRanges")
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) ‘GenomicRanges’

There is a binary version available but the source version is later:
binary source needs_compilation
GenomicRanges 1.30.0 1.30.3 TRUE

Binaries will be installed
trying URL 'https://bioconductor.org/packages/3.6/bioc/bin/windows/contrib/3.4/GenomicRanges_1.30.0.zip'
Content type 'application/zip' length 3223041 bytes (3.1 MB)
downloaded 3.1 MB

package ‘GenomicRanges’ successfully unpacked and MD5 sums checked
Warning: cannot remove prior installation of package ‘GenomicRanges’

The downloaded binary packages are in
Z:\Users\Marjory\AppData\Local\Temp\Rtmp8CRT7c\downloaded_packages

Then when attempting to install GenomicFeatures:

> ## try http:// if https:// URLs are not supported
> source("https://bioconductor.org/biocLite.R")
Bioconductor version 3.6 (BiocInstaller 1.28.0), ?biocLite for help
> biocLite("GenomicFeatures")
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) ‘GenomicFeatures’
also installing the dependency ‘GenomicRanges’

There is a binary version available but the source version is later:
binary source needs_compilation
GenomicRanges 1.30.0 1.30.3 TRUE

Binaries will be installed
trying URL 'https://bioconductor.org/packages/3.6/bioc/bin/windows/contrib/3.4/GenomicRanges_1.30.0.zip'
Content type 'application/zip' length 3223041 bytes (3.1 MB)
downloaded 3.1 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/bin/windows/contrib/3.4/GenomicFeatures_1.30.3.zip'
Content type 'application/zip' length 2737145 bytes (2.6 MB)
downloaded 2.6 MB

package ‘GenomicRanges’ successfully unpacked and MD5 sums checked
Warning: cannot remove prior installation of package ‘GenomicRanges’
package ‘GenomicFeatures’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
Z:\Users\Marjory\AppData\Local\Temp\Rtmp8CRT7c\downloaded_packages

I have tried installing through both RStudio and the RGui, tried opening and closing, tried the clearing of the envrionment too, but it is still telling me that there is no package called GenomicRanges or GenomicFeatures etc.

Answer: Please take a look at the comment by Martin Morgan, here: https://support.bioconductor.org/p/77688/#77693

> I suggest that as administrator you start R and run the command
> .libPaths(). Then quit R, remove the version of R that you have
> installed, and remove the directories reported by .libPaths().
Question: Sorry if this may appear like a really dumb question, but I am new to RNAseq and feeling lost about some of the terminologies and calculations used in literature.

I understand what's the difference between the principle behind paired-end and single-end reads,
however I need to clarify the method to calculate number of reads in paired end vs. single end reads.

For single end reads, I was told that I could calculate the no. of reads by:
wc -l , divided by 4.

For paired end reads, because 2 files are generated (R1, R2), so is it correct that I calculate the number of reads by adding the read count derived from R1 and R2?

Thank you very much for your clarification.

Answer: In paired-end sequencing its often less confusing to talk about fragment numbers rather than the more ambiguous read numbers. So count the number of reads in R1. That's the number of fragments sequenced (it's also the number of reads in R2). If you really prefer the words "reads", you could use "we sequenced XXX read pairs", so it's clear what's meant.
Question: **What is the code for male and female in plink covariate files? Is it (male = 1, female = 0), or ('1' = male, '2' = female, '0' = unknown)?** Is the code for gender different in covariate and fam. files?

I'm wondering about the inconsistent statements in the the manual of plink 1.9, as showed below. Thank you in advance.

---
As mentioned in [the manual of plink 1.9 which discusses the coding for male/female of covariate files][1]: "By default, when at least one male and one female is present, sex (male = 1, female = 0) is automatically added as a covariate on X chromosome SNPs, and nowhere else."

As mentioned in [the manual of plink 1.9 which discusses the coding for male/female of .fam files][2]: Sex code ('1' = male, '2' = female, '0' = unknown).

[1]: https://www.cog-genomics.org/plink/1.9/assoc#linear
[2]: https://www.cog-genomics.org/plink/1.9/formats#fam

Answer: The expected coding in the .fam file is male='1', female='2'; this is then coded by plink 1.x as male=1, female=0 during the --linear/--logistic regression. In other words, you would get the same results (outside of chrX) with `--linear sex` as you would with `--linear` combined with a covariate file with a male=1, female=0 sex covariate added.

Yes, this is counterintuitive, so I got rid of this discrepancy in plink 2.0. From its [--glm documentation][1]: "Note that PLINK 2.0 encodes the .fam/.psam sex covariate as male = 1, female = 2, to match the actual numbers in the input file. This is a minor change from PLINK 1.x." So with PLINK 2.0, if you use male=1, female=2 coding in both file types, you don't have to worry about the sign of the sex beta coefficient changing on you with .fam vs. `--covar`.

(With that said, even with PLINK 1.x, you don't have to worry about any p-values being affected by the 1/2 vs. 1/0 coding.)

[1]: https://www.cog-genomics.org/plink/2.0/assoc#glm
Question: Hello everyone!

I use apache (v2.4.18) on a GNU/linux (ubuntu) machine. /var/www/ is configure as the root for the homepage, jbrowse is installed in /var/www/html/summitdb/jbrowse, however the data is too big to store there.

My question is: how can I reach data with jbrowse outside /var/www/?

I have already tried using an alias without success:

I added the following to /etc/apache2/sites-available/000-default.conf

Alias "/data3" "/data/projects/summitdb/jbrowse/data3/"

Require all granted

Copied the trackList.json inside /var/www/html/summitdb/jbrowse/data and modified it like so:

"urlTemplate" : "/data3/tracks/100206/{refseq}/trackData.json",

The error massage is this:
SyntaxError: missing : after property id when parsing fsl/trackList.json.

Actually this using an alias is recomanded at http://gmod.org/wiki/JBrowse_Configuration_Guide#Configuring_track_locations_with_Apache.
What am I missing?

I also tried simply creating a symbolic link to my data and reach that by giving th link as atribute to the url (www.myhompa.ge/summitdb/jbrowse/index.html?data=symlink) but that does not seem to work either. Even though apache2 is set to allow symlinks.

In this case the error message is this:
RequestError: Unable to load http://myhompa.ge/summitdb/jbrowse/symlink/trackList.json status: 403.

This is my first time with both with jbrowse and apache2 and I am very confused. Any help would be appreciated!

Thanks,
Levi

Answer: OK, it's working now. There were actually two tings. There was a line uncommented in the json file that caused syntax error (possibly my metadata adding script is not working perfectly). The other one was solved by moving the jbrowse out of /var/www/ to the folder where the data is, and making a symlink there. To me, it seams like apache2 actually follows the symlink but jbrowse not.

Thanks to anyone who gave it a thought!

Chears,
Levi
Question: Hi
I have many large files with association results. Each file containes 8 columns (the 3rd one is the p value), I need to create from each file a new one conataining only observations where the p value is < 10 e-5.
How can I do this with bash code ?
Here a small example from these files :

SNP N P p2 or1 or2 q q1
c10_pos5974849 2 0.1881 0.1881 1.1931 1.1931 0.5707 0.00
c10_pos5975482 2 0.3225 0.3225 0.8670 0.8670 0.8840 0.00
c11_pos68438345 2 0.6537 0.66 0.9705 0.9690 0.2856 12.29
c11_pos107693921 2 0.8938 0.8558 1.0133 1.0250 0.1755 45.52
c12_pos67499221 2 0.8351 0.8351 1.0236 1.0236 0.6413 0.00
c14_pos67844869 2 0.1103 0.1915 0.7334 0.7229 0.2039 38.05
c14_pos68073026 2 0.09954 0.1298 0.6383 0.6215 0.2662 19.11
c14_pos68087872 2 0.3704 0.3704 1.2500 1.2500 0.7319 0.00

Thank you

Answer: find . -type f -name "*.common.suffix" | while read F ; do awk '($5=="P" || $5 < 1E-5)' $F > ${F}.subset.txt ; done
Question: Hi Biostars Users,

I am rookie in bio-informatics interested in undergoing NGS training to learn how to build NGS pipeline for whole exome sequencing analyses.

May you please guide me to find a suitable training program?

Best,
Kizuna

Answer: Perhaps one of out trainings might be of interest to you!

Check out http://www.ecseq.com/training.html
Question: Hello,

I am currently focusing on identifying denovo mutations from my trio data (parents are unaffected and child is affected). I used PhaseByTransmission. However, I found all denovo mutations (child is heterozygous, and both parents are hom. ref) were not phased (i.e. I am getting '/' instead of '|'). Do you think it is an error? If I search autosomal recessive, variants were phased correctly. What is the problem in my analysis? I am pasting the summary results provided by PhaseByTransmission for your kind perusal. Please also comment on the summary results, are they looking odd?

Please help.

```
java -jar /gatk_3.3/GenomeAnalysisTK.jar -R /reference_sequence/human_g1k_v37.fasta -T PhaseByTransmission -V trio1.vcf -ped trio1.ped --DeNovoPrior 0.00001 -o trio_out.vcf --MendelianViolationsFile mendelian_violation.vcf

INFO 20:04:04,201 GenomeAnalysisEngine - Strictness is SILENT
INFO 20:04:04,341 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000
INFO 20:04:04,453 PedReader - Reading PED file trio1.ped with missing fields: []
INFO 20:04:04,457 PedReader - Phenotype is other? false
INFO 20:04:04,510 GenomeAnalysisEngine - Preparing for traversal
INFO 20:04:04,530 GenomeAnalysisEngine - Done preparing for traversal
INFO 20:04:04,531 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
INFO 20:04:04,531 ProgressMeter - | processed | time | per 1M | | total | remaining
INFO 20:04:04,532 ProgressMeter - Location | sites | elapsed | sites | completed | runtime | runtime
INFO 20:04:34,824 ProgressMeter - 15:96876611 147844.0 30.0 s 3.4 m 77.5% 38.0 s 8.0 s
INFO 20:04:43,701 PhaseByTransmission - Number of complete trio-genotypes: 139299
INFO 20:04:43,702 PhaseByTransmission - Number of trio-genotypes containing no call(s): 0
INFO 20:04:43,703 PhaseByTransmission - Number of trio-genotypes phased: 124651
INFO 20:04:43,703 PhaseByTransmission - Number of resulting Het/Het/Het trios: 13391
INFO 20:04:43,704 PhaseByTransmission - Number of remaining single mendelian violations in trios: 937
INFO 20:04:43,704 PhaseByTransmission - Number of remaining double mendelian violations in trios: 12
INFO 20:04:43,704 PhaseByTransmission - Number of complete pair-genotypes: 0
INFO 20:04:43,705 PhaseByTransmission - Number of pair-genotypes containing no call(s): 0
INFO 20:04:43,705 PhaseByTransmission - Number of pair-genotypes phased: 0
INFO 20:04:43,705 PhaseByTransmission - Number of resulting Het/Het pairs: 0
INFO 20:04:43,706 PhaseByTransmission - Number of remaining mendelian violations in pairs: 0
INFO 20:04:43,706 PhaseByTransmission - Number of genotypes updated: 4395
INFO 20:04:45,481 ProgressMeter - done 201351.0 40.0 s 3.4 m 100.0% 40.0 s 0.0 s
INFO 20:04:45,482 ProgressMeter - Total runtime 40.95 secs, 0.68 min, 0.01 hours
INFO 20:04:47,002 GATKRunReport - Uploaded run statistics report to AWS S3
```

Answer: Those are your likely candidate de novo mutations, the tool uses `/` instead of `|` to indicate a violation of mendelian inheritance pattern and the presence of an un-inherited allele at these loci.
Question: I was wondering which ways are available for getting metadata from accession numbers. I have seen other tools such as Nextstrain make use of a so called "metadata" file to describe used sequences. The file looks something like this:
![Metadata for sequences][1]

[1]: https://imgur.com/a/uL3m7T5
https://imgur.com/a/uL3m7T5

It shows various data from NCBI for the accession numbers such as virus strain, country, date, URL, etc.
For me the most import ones are strain, country and date. Are there ways to download such data automatically when you have a list of accession numbers?

Any help is appreciated.

Answer: Using [EntrezDirect][1] :

$ esearch -db nuccore -query "KY317939" | esummary | xtract -pattern DocumentSummary -element SubName
ZIKV/Homo_sapiens/Colombia/2016/ZC204Se|Homo sapiens|Colombia|serum|06-Jan-2016|Antibody Systems Inc

Fields you are getting above are (separated by `|`)

isolate|host|country|isolation_source|collection_date|collected_by

[1]: http://bit.ly/entrez-direct
Question: Dear Friends,

I mostly script in shell/bash as am comfortable with it and also use a lot of linux commands and Entrez utilities for various big file analysis; however, when I stared working on VCF files, dealing with million rows and columns I find issues in operating these files with bash. For instance, below is an example script to count number of SNPs for samples(TCGA...barcode...) in a VCF file:

for i in {10..10000}
do
snp_count=$(zcat xxx.vcf.gz | grep -v "##" | awk "{print \$$i}" | grep -v "TCGA" | grep "0/1\|1/1" | wc -l)
sample=$(zcat xxx.vcf.gz | grep -v "##" | awk "{print \$$i}" | grep "TCGA")
echo -e "$sample\t$snp_count\n"
done

Output:
TCGA...barcode1 60123
TCGA...barcode2 45245
.
.

The script works but it has been running for 2 days now, and will take time. I would really appreciate if you could share/suggest on the best ways to deal with VCF files using bash/shell or if another programming language could make a difference? However, if possible, I would prefer bash scripting ways to tackle VCF files.

Thanks,
DK

Answer: Here is oneliner using bioalcidaejdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

$ wget -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chrY.phase3_integrated_v2a.20130502.genotypes.vcf.gz" |\
gunzip -c |\
java -jar dist/bioalcidaejdk.jar -e 'stream().flatMap(V->V.getGenotypes().stream()).filter(G->!(G.isHomRef()||G.isNoCall())).map(G->G.getSampleName()).collect(Collectors.groupingBy(Function.identity(),Collectors.counting())).forEach((k,v)->println(k+" "+v));' -F VCF
(...)
NA20852 441
HG02489 1233
HG02008 112
NA11919 905
HG02002 119
HG01398 133
HG03577 837
HG02009 1217
NA18549 819

detail:

* `stream().` get the stream of variants
* `flatMap(V->V.getGenotypes().stream()).` convert it to a stream of genotypes
* `filter(G->!(G.isHomRef()||G.isNoCall())).` remove the hom_ref and the nocall
* `map(G->G.getSampleName()).` extract the sample name from the genotype
* `collect(Collectors.groupingBy(Function.identity(),Collectors.counting())).` count the number of distinct names
* `forEach((k,v)->println(k+" "+v));` print it
Question: I'm having a plink file (ped map format) which I have subset for a list of SNPs and now I want to know with which genes these SNPs are related.

I did as the website suggested,

plink --file mydata --assoc

so the `plink.assoc` file was created as should be, but then when I run

plink --noweb --gene-report plink.assoc --gene-list glist-hg18

the `plink.range.report` file which should contain what I want, is empty.

I have manually checked and confirmed that for at least one SNP there is one corresponding gene.

When I run

plink --noweb --gene-report plink.assoc --gene-list glist-hg18 --gene-report-empty

which is the opposite command, it contains all the genes of the glist-hg18 file.

At this point I am unable to think of a solution.

Answer: Can you check whether the latest PLINK 1.90 build has the same problem? (Note that it does not support the undocumented `--gene-report-empty` flag.) (Edit: oops, this should not have been an answer, but hopefully it will lead to one in short order.)
Question: Hello,

I was wondering what is the quickest way to get a listing of the human Gene Symbols from Refseq Build 37.
Thannks in advance for your suggestions.

Fred

Answer: curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz" |\
gunzip -c | cut -d ' ' -f 13 |\
sort -u
Question: Hello there,

I did standard GSEA for my microarray expression data. As I found out GSEA tools provide GSEA plot just for top 20 of enriched pathways. How can I provide plot for other pathways which are also significantly enriched (adj-pval < 0.05) in my dataset based on GSEA?

I am trying *fgsea* package in r. Could anyone please bring an example here that shows, how the input should look like for plotEnrichment()? What's wrong with my input?

I want to plot for TGFb gene set and this is my script : `plotEnrichment(tgfb_gene, gene_list, gseaParam = 1, ticksSize = 0.2)`

so for the first argument, *examplePathway* I used this list:

> tgfb_gene

$HALLMARK_TGF_BETA_SIGNALING

[1] "TGFBR1" "SMAD7" "TGFB1" "SMURF2" "SMURF1" "BMPR2" "SKIL"
[8] "SKI" "ACVR1" "PMEPA1" "NCOR2" "SERPINE1" "JUNB" "SMAD1"
[15] "SMAD6" "PPP1R15A" "TGIF1" "FURIN" "SMAD3" "FKBP1A" "MAP3K7"
[22] "BMPR1A" "CTNNB1" "HIPK2" "KLF10" "BMP2" "ENG" "APC"
[29] "PPM1A" "XIAP" "CDH1" "ID1" "LEFTY2" "CDKN1C" "TRIM33"
[36] "RAB31" "TJP1" "SLC20A1" "CDK9" "ID3" "NOG" "ARID4B"
[43] "IFNGR2" "ID2" "PPP1CA" "SPTBN1" "WWTR1" "BCAR3" "THBS1"
[50] "FNTA" "HDAC1" "UBE2D3" "LTBP2" "RHOA"

and for the second argument, *exampleRanks*, I used a **vector** of the **logFC** (output of differential gene expression analysis), the names of vector is the name of genes. But I got this error and warnings:

Error: Aesthetics must be either length 1 or the same as the data (1): x, xend
In addition: Warning messages:
1: In max(tops) : no non-missing arguments to max; returning -Inf
2: In min(bottoms) : no non-missing arguments to min; returning Inf
3: In max(tops) : no non-missing arguments to max; returning -Inf
4: In min(bottoms) : no non-missing arguments to min; returning Inf
5: In max(tops) : no non-missing arguments to max; returning -Inf
6: In min(bottoms) : no non-missing arguments to min; returning Inf

Is there any other ways to provide or re-plot GSEA results?

I really appreciate any help!!

Answer: I have been working on GSEA too and I found two very interesting resources but I am not that proficient with this kind of analysis so I am still not able to interpret the plots.

However, here you can find a [very nice tutorial for R][1] that gives you some plotting. Also there is a [question][2] and another [tutorial][3] about how to interpret the graphical results.

Finally, this is an [article][4] about GSEA.

Hope this helps!

[1]: https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html#gsea-analysis
[2]: https://www.biostars.org/p/367191/
[3]: http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.htmlhttp://
[4]: https://www.pnas.org/content/102/43/15545
Question: Hello Biostars Community,

General question(s) here, specifically in regards to sequencing platforms and some related questions on sequencing, as well?

1. Could using a newer/the newest computational genome/annotation (for example, presently, Ensembl 107 or the newest Gencode version) adversely effect the actual truth of what was sequenced?

2. When sequencing is done through an Illumina machine or other big name company machines, are those sequencing platforms *completely independent* from the genome or DNA/cDNA being sequencing?

3. What happens if, for example, "famous gene ABC" and "low-profile gene XYZ" are found to have different 3' and/or 5' ends by some new discovery, in Illumina, would adapters still link to them to perform those bridge PCR reactions on the flow cell lanes, or would it be that all the data published before on "famous gene ABC" and "low-profile gene XYZ" should be revisited? Or is it like question #2, "completely independent" - are even adapters independent from genes?

I was doing some reading, and I guess the gene sequence really only matters for probe-based sequencing (chips and arrays?). Hopefully this question could be a good resource for others?

Thank you in advance.
- Pratik

Answer: A newer genome annotation/curation would be **independent** and not affect the actually nucleic acids sequenced by the machine (unless they were captured using a targeted approach **prior** to sequencing - e.g. whole exome sequencing). Whole genome sequencing will include ALL fragments present (regardless of genome/annotation) - assuming they are physically "able to be sequenced" (see homopolymers, etc... in the comment above).

**Adapters are independent of genes**

> was doing some reading, and I guess the gene sequence really only
> matters for probe-based sequencing (chips and arrays?). Hopefully this
> question could be a good resource for others?

You are correct
Question: Hi,

I am performing a shuffle of reads from a ChIP-Seq experiment over some features in the genome and since I do not know the original distributions I am trying to perform an observed to expected calculation to see which features are enriched. Boiled down this means I have counts of reads mapping to feature x, then I perform a shuffle of the reads over the entire genome, count those and then calculate the count of how many reads mapped to feature x by chance. Unfortunately, I am not quite sure on how to get a p-value from this.

I hope you can help me with this.
Thank you very much!

Answer: This looks like a permutation test. Your p-value would be the number of times the permutations give more reads mapping to the feature under consideration divided by the number of permutations.
And if you're doing multiple tests (i.e. testing many features), you should do a correction for multiple testing. Bonferroni is fine but too conservative, you should prefer FDR.
Question: Hi all,
Hope someone can help.

We can pass a region to samtools view such as 1:1000-1010.

Is it possible to pass a list of regions to samtools view and have the output into a single sam file i.e. (something like)

$ samtools view -h .bam regionlist.txt > mysam.sam

what would the syntax for this be?

also, in what format does the region list file should be in?

thanks,
a

Answer: Samtools have the `-L` option

samtools view -b -L ROI.bed file.bam > ROI_file.bam

BUT `-L` does not use the samtools index so the search is slooow depending on how large is your bam file. In my benchmarks querying for 10, 100 or 1000 sequences take exactly the same time so seem that bed size does not matter too much. You need to do a bit of benchmarking if you prefer to do a `-L` or to do a loop (for a bam with 23 million reads take the same time the `-L` that writting all the stuff for looping ;-)) YMMV

$ time samtools view -b -L 100_ROI.bed file.bam > ROI.bam
real 0m38.831s
user 0m37.666s
sys 0m0.556s

$ time (samtools view -H file.bam > roi_xargs.sam; \
cat 100_ROI.bed | perl -lane 'print "$F[0]:$F[1]-$F[2]"' | xargs -n1 -t -I{} samtools view file.bam {} >> roi_xargs.sam; \
samtools view -bSh roi_xargs.sam > roi_xargs.bam \
)
real 0m7.188s
user 0m5.080s
sys 0m1.304s

And if you feel adventurous and your disks are using isilon or lustre you can use gnu-parallel and forget about same file concurrency and do it in parallel in a breeze. The only issue is that you would need to do the bam merge afterwards.