cannin/biostars_qa huggingface dataset:
Question: Hello,

I have a BED file that has many intervals that are overlapping.

My objective is to merge the intervals to see how much of the chromosome the bed file spans.

BEDtools merge was my natural go-to method, however, there's a catch.

I only want to merge overlaps that reciprocally overlap by say, 90%. BEDtools merge as far as I know, doesn't have this feature. Does anyone know any tool that I can use that can do this?

Thanks a bunch and merry Christmas!

Answer: With [BEDOPS][1] *[bedmap][2]*:

$ bedmap --count --echo-map-range --fraction-both 0.9 --delim '\t' intervals.bed \
| awk '$1>1' - \
| cut -f2- - \
| sort-bed - \
| uniq - \
> answer.bed

If the intervals are unsorted or are of unknown sort state, first [sort][3] before mapping:

$ sort-bed intervals.unsorted.bed > intervals.bed

On further thought, I think you may want to filter single-overlap cases; please see discussion further in this thread.

**Edit:** Changed `--fraction-either` to `--fraction-both` to do a correct mutual overlap test.

[1]: http://bedops.readthedocs.io/en/latest/index.html
[2]: http://bedops.readthedocs.org/en/latest/content/reference/statistics/bedmap.html
[3]: http://bedops.readthedocs.org/en/latest/content/reference/file-management/sorting/sort-bed.html
Question: AFIK, Short sequencing and short pairwise sequencing generated by next-generation-sequencing is common now-a-days and there exists several alignment tools for them.

But I want to ask about long pairwise sequencing. There exists some alignment tools for long sequences like [rHAT][1]. That means there is at least one tool which generates long reads.

Now, I am asking that is there any tool which generates **Long Pairwise Reads**? if yes, then is there exist any aligner or mapper for that? Give me some idea about it.

Thanks in advance.

[1]: https://github.com/derekguan/rHAT

Answer: See this nice graphic for various sequencing technologies currently available:
![enter image description here][1]

from: https://flxlexblog.wordpress.com/2016/07/08/developments-in-high-throughput-sequencing-july-2016-edition/

None of the long read technologies are paired-end. However, with long reads, the paired-end option becomes a lot less useful.

[1]: https://flxlexblog.files.wordpress.com/2016/07/developments_in_high_throughput_sequencing.jpg
Question: Hi everyone,

Sometimes I am asked to mention the most exciting moment or part of my research works, I often say that I was excited first time I know how to work with command lines in Linux. However I think might be there is a specific purpose behind asking this question. What would you reply if you asked so?

Thank you very much

Answer: Seeing my work in a published paper, be it figures or analysis or, ultimately, a citation to a paper where I was one of the authors. Seeing my first high-profile paper get used in a Coursera genomics class that I had signed up for was a genuine thrill.
Question: Hi,

why we performed cufflinks and what does the assembled transcripts mean? I read the protocol, if I am right that cufflink will provide a reference transcriptome on which the reads can be aligned and enables use to calculate the expression level of genes???

Answer: Cufflinks builds the transcriptome for a sample. Transcripts are assembled using aligned reads which you provide in the form of a bam file. It also quantifies the assembles transcripts. Reference transcriptome is provided if it is available for that species. It helps cufflinks in building the transcriptome of the sample.
Question: Hi Biostars!

I m meddling with one doubt.. Could anyone tel me why would we do loop modeling of nucleic acids?

Answer: Loop modeling in nucleic acids (i.e. in RNA) is done for the same reasons as loop modeling in proteins: to find out the possible 3D structure of certain regions in a biomolecule (here, RNA). For instance, a structure resulting from an X-ray experiment might miss a couple of bases/nucleotides, possibly because they couldn't be traced properly during processing (one example would be the PDB entry 1u9s of an RNAse P domain that misses the atomic coordinates for a complete tetraloop hairpin). Another use-case would be to check which loops could fit into a specified site in an RNA molecule, e.g. to find alternative loop structures that might serve (at least from a geometrical point of view) as replacement for an existing structure.

In any case, "loop modeling" in RNA is a 3D structure problem. Here, secondary structure prediction can only help you with finding the stem and loop regions in a sequence with unknown 3D structure. Since you already know the 3D structure of the molecule (with the exception of your target loop), secondary structure prediction does not yield any useful information for the loop modeling at all.

References:

Schudoma et al., 2010, Nucleic Acids Research
http://nar.oxfordjournals.org/content/38/3/970

Schudoma et al., 2010, Bioinformatics
http://bioinformatics.oxfordjournals.org/content/26/13/1671
Question: I am trying to install CopywriteR, which failed at updating a dependent package RMySQL. The error I got was as follows:

------------------------- ANTICONF ERROR ---------------------------

Configuration failed because libmysqlclient was not found. Try installing:

* deb: libmariadb-client-lgpl-dev (Debian, Ubuntu 16.04)

libmariadbclient-dev (Ubuntu 14.04)

* rpm: mariadb-devel | mysql-devel (Fedora, CentOS, RHEL)

* csw: mysql56_dev (Solaris)

* brew: mariadb-connector-c (OSX)

If libmysqlclient is already installed, check that 'pkg-config' is in your

PATH and PKG_CONFIG_PATH contains a libmysqlclient.pc file. If pkg-config

is unavailable you can set INCLUDE_DIR and LIB_DIR manually via:

R CMD INSTALL --configure-vars='INCLUDE_DIR=... LIB_DIR=...'

--------------------------------------------------------------------

My question is: where can I specify --configure-vars='INCLUDE_DIR=... LIB_DIR=...'? I cannot do this at command line because RMySQL is installed through bioconductor. If I do

R CMD INSTALL --configure-vars='INCLUDE_DIR=/mydir/mariadb-connector-c' RMySQL

I will get:

--------------------------------------------------------

Warning: invalid package ‘RMySQL’

Error: ERROR: no packages specified

--------------------------------------------------------

Thanks in advance!

Answer: I was able to install CopywriteR after the following steps:

1. Download CopyhelpeR_1.0.2.tar.gz

2. sudo R CMD INSTALL CopyhelpeR*.tar.gz

3. sudo R CMD INSTALL CopywriteR*.tar.gz

This somehow bypassed the RMySQL installation. Part of the problem might come from the R version: I was using R 3.2.3, while CopywriteR was built under R version 3.4.2.
Question: Hi guys, I have a simple question. I have RNA-Seq data from different batches. As suggested looking at many posts on-line I have pre-normalized my data (using the TMM from edgeR) then I have corrected them using Combat and then I have re-normalized them (for the library-size) using DESeq2. My question is: is it correct the second normalisation after Combat? Or at least is it not dramatically not-correct?

Thank you in advance

Best

Answer: A quick note since this is a common problem, but for batch correction you generally need to have multiple conditions per batch. If all your WT samples are in one batch, and all your KD samples are in another batch, you can't correct for it (as an example).

With that being said, you can usually add batch as a covariate to the regression formula in edgeR and DESeq2 as the simpler and more robust option. Your study design would look like the following example for DESeq2:

> df
condition batch
WT-1 WT batch_1
WT-2 WT batch_2
WT-3 WT batch_2
KO-1 KO batch_1
KO-2 KO batch_2
KO-3 KO batch_2

Your regression formula would then be `~ condition + batch`, which means your differential expression results for condition will be corrected for batch.
Question: Hi everyone,
I would like to know if there is any statistical motivation to plot -log(p-value) vs log(FC); why not -log(FDR) vs log(FC) directly?!?

Thanks! !!

Answer: There is none as far as I know reading different papers and reports, you can actually see [here][1] which reports volcano plot with -log10(FDR) vs log2(FC), at the end of the day it depends upon the user to use the best scaling metrics for the most represented visualization.

According to the definition it is the plots significance versus fold-change on the y- and x-axes, respectively. So your significance can be FDR corrected as well and in that case you are just restricting your FC values to a much stricter subsets to give more confidence to the visualization. There should not be any other statistical motivation other than giving it a more reliable plotting with less error prone significant points.

[1]: http://www.nature.com/articles/srep05698/figures/5
Question: Hello,

Is it possible to use vcftools weir-fst-pop on haploid genomes? I do not seem to be able to find any option to specify the ploidy of my samples.

Thank you in advance.

Answer: Fst involves measuring expected heterozygosity. Given that, it is impossible to calculate Fst for haploid organisms.
Question: Hello.

How to extract a protein sequence?

I have two files.

< file 1. complete_protein.fasta >

>protein_1

>DCXSTEISLFHEIWLF

>protein_2

>AJFOWIDJLSIDJFJ

>protein_3

>DJFLWIDJFLSKDJFL

>protein_4

>DKSJFLEISJDKJF

< file 2. only proteinID.fasta >

>protein_1

>protein_4

I need about sequence in file2. That sequence have in file1.

So, I tried "diff" command, but result is I don't want data.

How to extract??

Answer: Hi,

If you are working with linux distribution you can try the following command-line:

grep -f target_protein.txt -A1 protein_file.txt | sed '/--/d' > retrieved_protein.txt

This takes the `protein_file.txt` file that contains the same content that you posted above:

>
protein_1
DCXSTEISLFHEIWLF
protein_2
AJFOWIDJLSIDJFJ
protein_3
DJFLWIDJFLSKDJFL
protein_4
DKSJFLEISJDKJF

And the file `target_protein.txt` that contains the target protein names that you want to retrieve from the file above:

>
protein_1
protein_4

The output of protein target sequences of interest is saved in `retrieved_protein.txt`, that looks like:

>
protein_1
DCXSTEISLFHEIWLF
protein_4
DKSJFLEISJDKJF

I hope this helps,

António
Question: 1. How did it come to be that the alternate nucleotide was more frequent than the reference nucleotide?
2. How does one account for this phenomenon when designing a strategy to filter for variants of interest? Should I go through the complicated process of selecting those individuals who DO NOT have the variant and calculate that the REFERENCE frequency in the population is probably around (1 - esp6500siv_all)?

I am researching a rare disease and have whole exome sequence data with the corresponding variant calls. Each variant call has been passed to annovar and among other data, we have looked up the frequency of the variant in the esp6500siv2_all data. Clearly a variant that was observed to have a high frequency in our sample but that had low frequency in esp6500siv2_all would be of disproportionate interest.

Low and behold I was surprised to find that 13% of the all of our variants (4055 out of 32131) had an allele frequency that was greater than 0.5. How can that be? I expected that all the allele frequencies would be <0.5.

I had thought that the variants would be akin to a minor allele frequency (MAF). Clearly I was wrong. I pulled 3 random variants from among the variants that had more than 0.5 frequency, to check them against the Exome Variant Server.

avsnp147 Chr Start End Ref Alt Gene.refGene esp6500siv2_all
1: rs3803530 15 89632842 89632842 C A KIF7 0.5373
2: rs621383 3 125118840 125118840 T C SLC12A8 0.9988
3: rs633561 11 64229857 64229857 A G NUDT22 0.9418

Looking up at NHLBI Exome Sequencing Project (ESP) [Exome Variant Server][2] and using All Allele

1. rs3803530: C>A; A=6984/C=6014 which means A is 6984/(6984+6014) or 0.537
2. rs621383: T>C; C=12479/T=15 which means C is 12479/(12479+15) or 0.999
3. rs633561: A>G; G=12240/A=756 which means G is 12240/(12240+756) or 0.941

[1]: https://drive.google.com/file/d/1jA7Azm2f1ss_rw3NoVgMSJ75tBqAf_qt/view?usp=sharing
[2]: http://evs.gs.washington.edu/EVS/

Answer: reference ≠ major ≠ ancestral ≠ wildtype

As Kevin says, the reference is just whatever is in the reference sequence, which is the sequence of whoever they happened to sequence for that region.

GRCh38 is an improvement compared to GRCh37 because the GRC sought out some loci where the reference allele was not the major allele in the 1000 Genomes project, and replaced those regions with tiny contigs which did have the major allele. Some, not all.
Question: Hi !

I have many .bam that I want to get their .bai using samtools in the terminal.

I tried the following command :

samtools index *.bam

However, I did not get any .bai file.

Regards

Answer: Samtools index only accepts a single input file, so using a shell metacharacter to specify multiple files will not work. I usually use a shell wrapper to run samtools index on a single file at a time.

```
#!/usr/bin/env bash

# index_all.sh

for INFILE in "$@"
do
samtools index $INFILE
done
```

Then it is simple to run:

./index_all.sh /directory/*.bam
Question: Hey all,
searching for ideas and advices.
i have 4 of enrichment results for biological processes. They should be relatively different from each other but i expect a certain number of categories to overlap between each other. My idea was indeed plotting all results in a comprehensive way, in order to show the overlapping categories (storing their enrichment folds), and also showing the "private" categories from each result.
any idea or suggestion will be much appreciated

Answer: Something like this should be just fine (this is recent code that I pieced together for a project):

This is built on top of the *ComplexHeatmap* package in R, with much customisation by myself. If it looks good, then I can guide you further with code examples. I have masked the gene names just to protect information.

Kevin
Question: Hi everyone,

I am working with genome wide data (SNPs) using plink and want to compute pairwize LD (r²) between SNPs for each chromosome (29 chromosomes).

This is an example of script I used:

plink --bfile myData1 --r2 --ld-window-r2 0 --out chr1
plink --bfile myData2 --r2 --ld-window-r2 0 --out chr2
. . . . . .
. . . . . .
. . . . . .
plink --bfile myData29 --r2 --ld-window-r2 0 --out chr29

And the result is like :

CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2
23 10121 T2500005 23 14911 T2300007 1
23 10121 T2500005 23 15894 T2300003 0.0175439
23 10121 T2500005 23 43444 T2300006 0.416667
23 10121 T2500005 23 60163 T0015398 0.416667

In fact, it works for all chromosomes except chromosome 24, which give an embty table :

CHR_A BP_A SNP_A CHR_B BP_B SNP_B R2

The same probleme for 3 breeds (always no results for chr24).
Do you have an idea?

Tanks for help.

Answer: —chr-set *is* required for chromosome 24 here. Otherwise, it is treated as chrY.

It may also affect the results for chromosome 23. To be safe, you should always use it with nonhuman species.
Question: Hi,

Is there a pipeline/pacakage that can easily get me genotypes (eg in AA, AG, GG) type format from Affy genome wide SNP CEL files. Thus far I've played around with Birdseed and CLRMM. Both suffer from a combination of being impossible to install, having awful documention or producing unusable output. I literally want, CEL files go in, (annotated!!!!) genotypes and confidence scores come out. Any suggestions?

Answer: **Edit by Kevin Blighe on March 10, 2020**: Affymetrix was subsequently acquired by ThermoFisher, and the Genotyping Console can now be downloaded from HERE.

----------------
---------

Affymetrix genotyping console

For SNP 6 it supports

calling genotypes

copy number/LOH

copy segments data

copy number variation analysis

Also see manual
Question: Hi,

I am only interested in aligning DNAseq reads to certain genes. If I split my reference genome based on the coordinates of my gene of interest (as present in the GTF/GFF file) and then use BWA for aligning my reads to the resulting 'smaller reference genomes', will it be a good idea?

If yes, is there a threshold to the number of bases upstream and downstream of the gene coordinates that should be considered? And what caveats does this method involving splitting the reference genome can have that I should pay attention to?

My motto for using this method is to reduce alignment time as I am only interested in say 20-30 genes and not all genes.

Answer: > for aligning my reads to the resulting 'smaller reference genomes', will it be a good idea?

NO, you'll get false positives. It's the same as : https://www.biostars.org/p/4405/ (you're 'masking' a whole chromosome) . Citing Heng Li:

> This will lead to wrongly mapped sequences, spurious SNPs/indels calls and all sorts of problems. I cannot think of a single use case when masking [before mapping] may lead to better outcomes."
Question: This question was inspired by Peter Cock 's tweet

https://twitter.com/pjacock/status/118012750105546752

Most softwares (C, java...) use a int32 (unsigned or not) to store the length of the chromosome. It isn't enough when the length of a chromosome is greater than INT_MAX or UINT_MAX

# define INT_MAX 2147483647
# define UINT_MAX 4294967295U

So my question is:

Is there any resource where one can find the length of the chromosomes. Something like BioNumbers.

What's the length of the longest chromosome ?

Answer: With the publication of [A chromosome-based draft sequence of the hexaploid bread wheat (*Triticum aestivum*) genome][1], the initial release of the Wheat ta 3B chromosome is 774 Mbp (or to be exact, 774434471bp in file ``https://urgi.versailles.inra.fr/download/wheat/3B/ta3bPseudomolecule.genom.fa.gz``)

[1]: http://dx.doi.org/10.1126/science.1251788
Question: My task is to repeat the DATA analysis of RNA-seq data as presented in a journal article using the tophat cufflinks pipeline.

For simplicity Ill just mention the 4 controls

The authors run cufflinks without a reference annotation on each control "to detect possible novel transcripts" --> then cuffmerge on the results --> they then say they run cufflinks again using the merged transctiprts.gtf as the reference annotation. It seems over complicated.

Cufflinks requires a .BAM file as input but cuffmerge output doesnt give a BAM file....so the only way i can see they did it is by re running cufflinks on every sample for a second time (waste of time?) except this time using the cuffmerge output as the reference annotation. This would mean re running cuffmerge again also afterward.

Surely " to detect possible novel transcripts" doesnt require running cufflinks on everything twice....I mean, isnt this the whole point of cufflinks.

Thanks in advance.
Kenneth

Answer: The first Cufflinks run is to generate a new annotation for each sample to discover novel transcripts.
The Cuffmerge run is to merge together all the annotations for each individual sample to create one merged annotation of better quality.
The second Cufflinks run is to quantify the transcripts based on the merged annotation file.

Yes, it is complicated, and the results will contain many false positives.
More importantly, it's generally a waste of time, unless you're working on a poorly annotated genome.
For well-annotated genomes like the mice, human, or drosophila genomes, you shouldn't bother trying to discover novel transcripts.
Just use the most recent annotation available.
Question: Hello everyone,

I'm trying to do a simple t-test on my microarray sample in R. My sample looks like this:

gene_id gene sample_1 value_1 sample_2 value_2
XLOC_000001 LOC425783 Renal 20.8152 Heart 14.0945
XLOC_000002 GOLGB1 Renal 10.488 Heart 8.89434

So the t-test is between sample 1 and sample 2 and my code looks like this:

ttestfun = function(x) t.test(x[4], x[6])$p.value
p.value = apply(expression_data, 1, ttestfun)

It gives me the following error: Error in t.test.default(x[6], x[8]) : not enough 'x' observations
In addition: Warning message:
In mean.default(x) : argument is not numeric or logical: returning NA

What am I doing wrong? Please help.

Many thanks.

Answer: I think there's some misconceptions operating here from the original questioner. First and foremost, a t-test is not just a way of calculating p-values, it is a statistical test to determine whether two populations have varying means. The p-value that results from the test is a useful indicator for whether or not to support your null hypothesis (that the two populations have the same mean), but is not the purpose of the test.

In order to carry out a t-test between two populations, you need to know two things about those populations: 1) the mean of the observations and 2) the variance about that mean. The single value you have for each population could be a proxy for the mean (although it is a particularly bad one - see below), but there is no way that you can know the variance from only one observation. This is why replicates are required for microarray analysis, not a nice optional extra.

The reason a single observation on a single microarray is a bad proxy for the population mean is because you have no way of knowing whether the individual tested is typical for the population concerned. Assuming the expression of a given gene is normally distributed among your population (and this is an assumption that you have to make in order for the t-test to be a valid test anyway), your single individual could come from anywhere on the bell curve. Yes, it is most likely that the observation is somewhere near the mean (by definition, ~68% within 1 standard deviation, see the graph), but there is a significant chance that it could have come from either extreme.

Finally, I've read what you suggest about the hypergeometric test in relation to RNA-Seq data recently, but again the use of this test is based on a flawed assumption (that the variance of a gene between the 2 populations is equivalent to the population variance). Picking a random statistical test out of the bag, just because it is able to give you a p-value in your particular circumstance is almost universally bad practise. You need to be able to justify it in light of the assumptions you are making in order to apply the test.

BTW, your data does not look like it is in log2 scale (if it is, there's an ~32-fold difference between the renal and heart observations for the first gene above) - how have you got the data into R & normalised it?