cannin/biostars_qa huggingface dataset:
Question: Hi - when using VEP to get allele frequency information for a list of SNPs, I ask for the gnomad_AF column. I have noticed however that around a third of my list return an NA value. How should I intepret this? Is that of all combined frequency consoritums across multiple continents, this SNP has never been seen (i.e. Ultra Rare)? Some of the SNPs do appear in dbsnp however, but I guess the exomes/genomes used for gnomad might not coincide with those reported to dbSNP? This seems about right to me - any lab can report to dbSNP, but only subset of the population end up in consotriums such as 100k genomes and exac.

Could anyone provide a definitive answer on this? Can I count them as "ultra rare" if a SNP does not appear in gnomad?

Thanks

Answer: I think it does count as ultra rare if gnomAD doesn't have it. As long as the comparing tool is configured well in comparing multi-allelic sites (as in, you're *sure* the variant is not in gnomAD), it can be labelled ultra-rare.
Question: It might be a very stupid question for many of you but, since it's my first variant calling, I didn't figure it out yet.

I have **mpileup**'ped two bam files from two samples, then I filtered the results with **vcfutils.pl** and called the genotypes with **bcftools** **call**. Now I have a **VCF** file containing what I want, but I managed to analyze differences only between the two samples and the reference (which is an assembly coming from a line that is *different* *from both* *samples*).

I have the variants between my two samples and the assembly, but what if I want to detect the differences **between** the two samples? I did it with awk / sed / cut and other command line tools but is there maybe a better and more straight-forward way to do that? Perhaps using bcftools? Up to know, I didn't find it.

Any suggestion appreaciated!

EDIT:
I did it with bcftools gtcheck as well, that works fine. I am asking if there are **other** reliable tools to test!

EDIT 2:
As this post has many views now, I guess it will be useful for everyone to know that I used **bcftools isec** and it worked brilliantly.

Answer: I was recently in a very similar situation. I wrote a bash script that finds the "similarity" between two VCF files. Someone might find it useful:

#!/usr/bin/env bash

sort -u <(grep -v '^#' $1 | cut -f1,2,4,5) > a
sort -u <(grep -v '^#' $2 | cut -f1,2,4,5) > b
comm -23 a b > a_only
comm -13 a b > b_only
comm -12 a b > ab

numer=`cat a_only b_only | wc -l`
denom=`cat ab a_only b_only | wc -l`
dist=`echo "$numer/$denom" | bc -l`
sim=`echo "(1-$dist)*100" | bc -l`

echo $sim

# CAREFUL
rm -f a_only b_only a b ab

If you put it in `sim.sh` then you can use it like so:

`$ ./sim.sh sample_1.vcf sample_2.vcf`
Question: How to use step by step the software of IPA (Ingenuity Pathway Analysis) to analyze the interaction between these genes:

IL3RA
CD38
CS
MS4A1
GPC2

in the cases of **lymphoma** disease.

Thank you for your enlightenment.

Best regards,

Dito Anurogo
(A Ph.D. student at Taipei Medical University, Taiwan)

Answer: There are several youtube videos, manuals and videos on product owner's website. Please consult them. I think they might have training sessions by vendor. Most of these vendors employ application scientists and contact them.
Question: This is a duplicate of:

1. https://www.biostars.org/p/7372/ from 5.4 years ago and
2. https://www.biostars.org/p/65920/ from 3.5 years ago.

Nevertheless, I found the answers on the second post enlightening, and 3.5 years later it seems worthwhile to have an update from the community (and @lh3) on how you deal with these issues.

Notably, @Aaronquinlan's gemini from the second most upvoted answer to question #2 is still under active development. Are people using it? What else are people doing?

Answer: This is called "variant warehousing" and there are several open source and commercial efforts in various stages:

Golden Helix VSWarehouse
Paradigm4 SciDB
WuXi NextCODE
CMH Variant Warehouse
ViaGenetics Genesis
Curoverse Lightning
Intel GenomicsDB
Cloudera OMICS
Question: EDIT: SOLUTION

As @rpolicastro said, the problem had to do with the conda channels. The solution for me was:

conda create -n cutadapt -c conda-forge -c bioconda cutadapt python=3.9

ORIGINAL POST

Hello everyone

I'm having trouble installing cutadapt.

I just installed conda and created a new environment for cutapdat with a specific python version (since it send me a message before about that).

conda create -n cutadapt python=3.9

When I enter the command `conda install -c bioconda cutadapt` the following message appears

```
UnsatisfiableError: The following specifications were found to be incompatible with each other:

Output in format: Requested package -> Available versionsThe following specifications were found to be incompatible with your system:

- feature:/linux-64::__glibc==2.35=0
- feature:|@/linux-64::__glibc==2.35=0
- cutadapt -> libgcc-ng[version='>=9.3.0'] -> __glibc[version='>=2.17']

Your installed version is: 2.35
```

But when I run `conda list` in my cutadapt environment to see the packages it says:

```
packages in environment at /home/camila/miniconda3/envs/cutadapt:

Name Version Build Channel
_libgcc_mutex 0.1 main
_openmp_mutex 4.5 1_gnu
ca-certificates 2022.4.26 h06a4308_0
certifi 2021.10.8 py39h06a4308_2
ld_impl_linux-64 2.35.1 h7274673_9
libffi 3.3 he6710b0_2

libgcc-ng 9.3.0 h5101ec6_17
libgomp 9.3.0 h5101ec6_17
libstdcxx-ng 9.3.0 hd4cf53a_17
ncurses 6.3 h7f8727e_2
openssl 1.1.1o h7f8727e_0
pip 21.2.4 py39h06a4308_0
python 3.9.12 h12debd9_0
readline 8.1.2 h7f8727e_1
setuptools 61.2.0 py39h06a4308_0
sqlite 3.38.3 hc218d9a_0
tk 8.6.11 h1ccaba5_1
tzdata 2022a hda174b7_0
wheel 0.37.1 pyhd3eb1b0_0
xz 5.2.5 h7f8727e_1
zlib 1.2.12 h7f8727e_2
```

So I don't know what the problem may be. Any thoughts?

Thanks a lot.

Answer: It could be related to your channels or the solver is having a problem with your specific linux install.

First try ensuring proper channel order when you go to create the environment.

```
conda create -n cutadapt -c conda-forge -c bioconda cutadapt python=3.9
```

If that doesn't work try using the beta libmamba solver which is more powerful than the current one.

```
conda update -n base conda
conda install -n base conda-libmamba-solver

conda create -n cutadapt --experimental-solver=libmamba -c conda-forge -c bioconda cutadapt python=3.9
```
Question: Struggling with biostats and association study designs.

I initially wanted to do an association study comparing two populations and seeing which SNPs are significant. For example, I am looking at just variation between centenarians (people who live >100years) and a control group. Should I be including age as a covariate? I am interested in detecting longevity variants or anything that suggests a difference from the control group. I do not think age would be necessary since it is not a confounder because it is not an independent covariate. I think I was reading that adding independent covariates can decrease power of the study.

When do you guys start considering modeling with linear or logistic regressions instead of GWAS? Is it when you have a dependent variable and a predictor variable you are interested in? Would you guys think I should have added covariates?

Or is another way possible to change the case and control population to groups that reflect covariate status? For example, If I was interested in centenarians with Alzheimer's compared to a control population with Alzheimer's, a logistic regression with Alzheimer status be more appropriate than running a GWAS on them?

Sorry for so many questions.

Answer: I'm a bit self-taught in this but since no-one has answered in 2 days I'll give it a try, others can feel free to chip in.

>Should I be including age as a covariate?

You add covariates when you expect that they have an influence on your phenotype - for example, gender often has an influence, but most commonly in humans it's population structure (like PC1 and PC2 and sometimes PC3 from PCA software). In your case with old age I wouldn't use age as a covariate as it's practically identical to your phenotype of interest, it's not independent.

Generally, I'd add covariates if you have surprisingly tiny p-values, and especially if your QQ-plots look bad. For a good introduction to reading and interpreting QQ-plots, see here: Behavior of QQ-Plots and Genomic Control in Studies of Gene-Environment Interaction. How to plot QQ-plots depends on the format of your output (some software like GAPIT does it automatically).

Have a look at these two tutorials for more on covariates, QQ-plots and population stratification: http://sites.tufts.edu/cbi/files/2013/01/GWAS_Exercise4.pdf and http://sites.tufts.edu/cbi/files/2013/02/GWAS_Exercise6_Stratification.pdf (these also have code for plotting QQ-plots for PLINK results)

>When do you guys start considering modeling with linear or logistic regressions instead of GWAS?

I'm not sure I understand this question, regressions are a part of GWAS. In my experience running a regression with PLINK gives relatively similar results to a mixed linear model as for example implemented in GAPIT. The actual p-values are different but the SNPs with the lowest p-values stay roughly the same.
Question: I have an article in which two types of multiple protein structures docking are described (multiple docking and ensemble docking) but it is not clearly stated how these approaches differ.

Answer: Ensemble docking refers to the approach of generating multiple structures of the target molecule using molecular dynamics simulation and study docking to this ensemble of target structures. Multiple docking refers to the breaking down of the target structure into smaller components that are then used for docking in various combinations. Both of these approaches try to account for conformational changes of the target protein.
Question: I'm using the SRA toolkit to convert some SRA files to Fastq format. I've been looking at the documentation to make sure I'm doing things right, and the word spot keeps coming up. My question is twofold.

What is a spot and how does it differ from a read?

Where is this (officially) documented (or is it)?

The reason I've separated these two questions is that I think I know the answer to the first one, but I'm not sure and I can't find the answer in any of the documentation or online. Also, I expect more people will know the answer to question #1 than question #2.

Answer: This is the description I received from the SRA staff (Adam Stine).

The spot model is Illumina GA centric. The flowcells have the locations where the adapters have stuck them to the glass of the lane. There are X and Y coordinates that identify these 'spots'. As the camera reads the fluorescent flashes during sequencing, the coordinates indicate which spot the new base is added to. All of the bases for a single location constitute the spot. There may be one or more divisions of those bases for technical reads (adapters, primers, barcodes, etc) and there will always be at least one biological read (forward, reverse). I usually think of the technical reads as the "known" sequence and the biological as the "unknown". When we store the data, the bases for a single spot are all stored as one string with the description of where the breaks occur as well as the type of read each segment represents. The spot length is the expected total length for all reads (used as a check to make sure we have all the data). As an example, a 2x150 run with a 6bp barcode and 12bp primer on the forward read would have 4 reads.

0 - barcode basecoord 1

1 - primer basecoord 7

2 - forward basecoord 19

3 - reverse basecoord 151

---

But you only need to explain SRA about the barcode and primer is you submit sequences that contains it..In my case, a third party provided me with the BAM files and I do not have the untrimmed sequences.

So the SPOT datamodel is useful for supplying untrimmed BAM.. yet, enable you to specify where the biological reads begin.

In my case, I have 2X100 bp without index and I am only supplying the Application read with the adapter trimmed. so I simply submit.

0 - forward basecoord 1 (Application read)

1 - reverse basecoord 101 (Application read)

---
Question: I understand that different Illumina sequencers can output different numbers of reads (i.e. HiSeq instruments can produce 40-400 million reads/lane, whereas MiSeq can produce 5-25 million reads/lane). My question is why? What factors in the Illumina next-generation sequencing technology impact read count outputs from these sequencers? For example, can certain instrument settings increase or decrease the number of reads produced (assuming that equal amounts of DNA were supplied to the sequencer)?

Answer: The two factors that limit the amount of data are 1) the maximum cluster density (i.e., number of clones) that can be imaged and resolved by the camera, and 2) the surface area of the flow cell lane that's imaged (MiSeq flow cells are much smaller than HiSeq).
Question: I would like to use snakemake to analyze my data sets. As I am going to work with different organisms, I would like snakemake to create a folder for each of them when indexing the genome.

let's say I would like to work with human and mouse data.

In my `config.yaml` file I have the following snippet:

organism:
"Dmel"
# "Dpse"

Each time I am working with an organism I can comment out the others in the list.

Below are the rules for getting the fastA files as well as the annotation and the indexing step (as an example) with STAR.

#### get reference genomic data for the mapping (fasta and gtf files). Links are added in the config file
rule get_genome:
output:
fastA="genome/genome.fa",
gtf="genome/genome.gtf"
shell:
"""
wget -nc -O - {config[fastA]} | gunzip -c - > {output.fastA}
wget -nc -O - {config[gtf]} | gunzip -c - > {output.gtf}
"""

### Indexing the reference genome

rule star_index:
input:
fasta="genome/genome.fa",
gtf="genome/genome.gtf"
output:
directory("genome/starIndex/")
threads: 16
params:
prefix = {config["organism"]}
shell:
"mkdir -p {output} && "
"STAR --runThreadN {threads} "
"--outFileNamePrefix {output}{params.prefix} "
"--runMode genomeGenerate "
"--genomeDir {output} "
"--limitGenomeGenerateRAM {config[RAM]} "
"--genomeSAindexNbases {config[SAindex]} "
"--genomeFastaFiles {input.fasta} "
"--sjdbGTFfile {input.gtf} "
"--sjdbOverhang 100"

I would like to know how to add the value from the config file in the organism, e.g. `config[organism]` into the rules, so that, when I download the fastA and gtf they will be renamed accordingly.

I have tried to add `config[organism]` to the path, but it gives me an error. I have also tried this for example:

rule get_genome:
output:
fastA="genome/{config[organism]}.fa",
gtf="genome/{config[organism]}.gtf"

But it just creates the files `{config[organism]}.gtf` and `{config[organism]}.fa`.

Can someone please help me understand how I can add the parameter from the config file to the different paths in the various rules?

thanks

Answer: You'd need `wildcards`. https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#wildcards

[~/Data/scratch/tmp/biostar]$ cat wildcards.py
rule:
input: expand("{organisms}/starIndex/", organisms=config['organism'])

rule get_genome:
output: fa = touch('{organism}/{organism}.fa'),
gtf = touch('{organism}/{organism}.gtf')
run:
# remove `touch` and implement logic to download fa and gtf
# you can access organism via wildcards, i.e. config[wildcards.organism]
pass

rule star_index:
input: fa = '{organism}/{organism}.fa',
gtf = '{organism}/{organism}.gtf'
output: directory("{organism}/starIndex/")
shell: 'mkdir -p {output}'

[~/Data/scratch/tmp/biostar]$ snakemake -s wildcards.py --config organism=GRCh38
...
[~/Data/scratch/tmp/biostar]$ tree GRCh38/
GRCh38/
├── GRCh38.fa
├── GRCh38.gtf
└── starIndex

1 directory, 2 files

In this case, `config['organism']` can be a list of organisms and the workflow will expand itself. Hope the example helps.
Question: Are there any programs available (Python, Perl, whatever), that can take a degenerate nucleotide sequence and translate it into its multiple possible oligos?

Any help would be appreciated.

Answer: Like Pierre said, this is a combinatorial problem that can quickly blow up in your face. However, that aside, it easy to generate all the combinations using a recursion. Here is a compact implementation in Perl:

```
#!/usr/bin/perl -w

# Lookup table of degenerate IUPAC nucleotide codes.
my %deg2nuc = (
"R" => ["A", "G"],
"Y" => ["C", "T"],
"S" => ["G", "C"],
"W" => ["A", "T"],
"K" => ["G", "T"],
"M" => ["A", "C"],
"B" => ["C", "G", "T"],
"D" => ["A", "G", "T"],
"H" => ["A", "C", "T"],
"V" => ["A", "C", "G"],
"N" => ["A", "C", "G", "T"]
);

# Recursive function that replaces degenerate nucleotides with all combinations.
sub generate
{
if ($_[0] =~ /(.*)([RYSWKBDHVN])(.*)/) {
my $head = $1;
my $tail = $3;
my @seqs;
foreach my $nuc (@{$deg2nuc{$2}}) {
push @seqs, generate($head.$nuc.$tail);
}
return @seqs;
}
else {
return $_[0];
}
}

# Demo: print all sequences generated from ANCRG.
print join("\n", generate("ANCRG")), "\n";
```

EDIT:

On second thought, that was not at all compact by Perl standards. Here is the really compact version that ventures into *code golf* territory:

```
#!/usr/bin/perl -w

my %A = ("R"=>1, "W"=>1, "M"=>1, "D"=>1, "H"=>1, "V"=>1, "N"=>1);
my %C = ("Y"=>1, "S"=>1, "M"=>1, "B"=>1, "H"=>1, "V"=>1, "N"=>1);
my %G = ("R"=>1, "S"=>1, "K"=>1, "B"=>1, "D"=>1, "V"=>1, "N"=>1);
my %T = ("Y"=>1, "W"=>1, "K"=>1, "B"=>1, "D"=>1, "H"=>1, "N"=>1);

$_ = "ANCRG\n";
while (s/(.*)([RYSWKBDHVN])(.*)\n/(exists $A{$2} ? "$1A$3\n" : "").(exists $C{$2} ? "$1C$3\n" : "").(exists $G{$2} ? "$1G$3\n" : "").(exists $T{$2} ? "$1T$3\n" : "")/e) {}
print;
```
Question: Hi all,
My question is to find a method to count enhancer RNAs from RNASeq data. Although it may have been asked before, I could not find any post close to mine. I have a couple of fastq files of RNASeq data. I am using TopHat for alignment. When I put my .bam files to IGV, I can visually see some reads in locations where enhancers bind to. But when it comes to using a tool to count the reads in those locations, although those specific locations are included in my .gtf file with an Ensemble ID, the count is "0".

my code for featureCounts is
``../anaconda/bin/featureCounts -a .../bcbio_ku/share/bcbio/genomes/Hsapiens/GRCh37/rnaseq/ref-transcripts.gtf -o .../work/my_generated.counts -t transcript -s 0 -p -C .../final/Unif_DMSO/Unif_DMSO-ready.bam``

Is there anyone using an alternative method for enhancer RNA counting? I believe a lincRNA counter or miRNA counter would help me build the method too. The last option for me is to create a unique .gtf file according to "http://enhanceratlas.org/download.php" for each cell line through a python script.

Thanks already

Answer: You could use bedtools to count number of reads mapping to interest of your region

> intersectBed -wa -wb -s -a bedFile -b BAMFile

Check different parameters

Keep in mind that eRNA are generally not well detected in most RNA-seq samples because of the way how libraries are prepared.

1) As eRNA are generally non-polyadenylated, libraries should be prepared with ribo-depletion kit to make sure non polyadenylated RNAs are in the sample. If poly-A is used in library prep, eRNA are removed/depleted.

2) eRNA are degraded rapidly, so libraries should be sequenced at very high coverage to detect these events. Alternatively, instead of steady mRNA (like RNA-seq), sequencing the nascent RNA (like GRO-seq) enriches eRNA.

3) Not usually done, but if libraries are further size selected, lets says, 50-500 nt long, expected size range of eRNA, you will enhance the eRNA signals.
Question: Hello everyone,

I am interested in analyzing some data from TCGA, specifically the LAML project. From there I managed to download the HTSeq counts for the whole cohort, however I would like to map each sample to specific mutations for given set of genes in order to perform DE analysis between pools of samples with the same gene mutated.

I am struggling with this second part. I found a R package named RTCGAToolbox in which you can download the expression data and the Mutation data among others. Instead, what I am getting are the mutations present in the cohort rather than a table to map samples with their mutations. Is there any possibility to download this kind of metadata, or I have to manually check from MAF files for each sample? If I have to, how should I do it?

Thanks in advance to everyone,

Esteve

P.S: Checked before writing this post if there was any similar one, but I didn't find any. Apologies if it is a repeated question.

Answer: I did a lot of the analysis for the TCGA AML paper. The simplest way to pull mutations is to go to the [paper's supplemental site](https://tcga-data.nci.nih.gov/docs/publications/aml_2012/) and pull Supplemental Table 6. Note that the variants are all on build 36, so liftover may be in order if you're working from a newer reference genome build.
Question: Hello,

I've watched a video from YouTube about OMIM. Then I get that OMIM already had more than 3000 disease-associated genes till last year. Now I want to get all of those genes from OMIM to do some experiments.

I think maybe I should use "Gene Map" of OMIM, but I still don't know how to search to get all disease-associated genes.

Do you have any good advice or guide?

Thank you in advance,

naulty

Answer: I endorse the alternative BioMart. I'd perhaps point to another [link][1] as well. The reason for that is due to the fact that there can be a lag between the latest version of the Ensembl database from our main site (as part of our [release cycle][2]) and the version available in Central BioMart. At the moment they are in sync but as soon we release a new version, it takes some time for them to update it. So just keep an eye on that.

[1]: http://www.ensembl.org/biomart/martview/a718a2c2ac2c280d41fa088bc5595bd9
[2]: http://www.ensembl.info/blog/2014/04/17/ensembl-release-cycle-and-grch38/
Question: Does anyone know of any tools/scripts that would score a MSA in a sliding window manner?

I want to get the average identity across the length of an operon MSA, so I'm looking for something that gives me back a measure of conservation that I can graph.

I found some old threads on BioStars, but a lot of the linked packages either aren't exactly what I'm looking for, or the links etc are just defunct.

If I have to write this myself (probably with python), does anyone have any thoughts on appropriate ways to score each window?

Answer: Wrote my own in the end. Uses some bits of code I found online, but a bit more performant.

https://gist.github.com/jrjhealey/130d4efc6260dd76821edc8a41d45b6a
Question: My lab has sequenced an exome using the "Agilent SureSelect All Exon Kits" ( http://www.genomics.agilent.com/CollectionSubpage.aspx?PageType=Product&SubPageType=ProductDetail&PageID=3041 )

Does anyone knows where I can find a BED describing the captured regions ?

P.

Answer: According to this post http://www.biostars.org/p/5187 , maybe you can find the information in Agilent's Earray https://earray.chem.agilent.com/earray/
Question: HI,

Can someone give a good definition of the difference ?

Thanks

Answer: UCEs are short- and medium-sized stretches of DNA that are identical between different species. They are often non-coding.

Orthology refers always to coding sequences. It requires that function be conserved between different species, but there is no requirement for identical sequence.

While presumably there is a requirement for functional conservation in both cases, evolutionary pressure for sequence conservation is larger on UCEs than orthologs.