cannin/biostars_qa huggingface dataset:
Question: I routinely expect that a 1GB transfer of data across the internet to be done without a single bit flip.

As a bioinformatic newbie, I have no reference point for how repeatable genomic sequencing currently is. This generates a series of questions:

1) If I have my genome fully sequenced at two different, state of the art sequencing centers (e.g. the best in america, the best in europe) how similar would the results be (and how should we define and measure similarity)? (I'm interested in reference genome quality sequences, not e.g. clinical quality)

2) For that matter, if I use the same center, but just on two different occasions separated by a a couple of weeks, what should I expect?

3) Instead of human genome, suppose a strain of E. coli was being analyzed -- how does that effect the repeatability?

4) Naively I would assume that I would receive back data organized in 46 contiguous chunks. Pragmatically, what should I expect? -- e.g. you'll get haploid data, and the avg # of contigs per chromosome will be x.

5) Genomic Coordinate Systems: For me, similarity between two sequencing runs means first, that there is computational simple way to translate the coordinates of a region of interest in one data set to the other, and second, the sequence data matches. For instance, how does this work for humane genome GRCh37 vs GRCh38 -- how would I go about translating coordinates for e.g. all the opsins in GRCh37 to their coordinates in GRCh38, and how likely would it be that there would be a single nucleotide mismatch for this limited portion of dna?

Hope those questions make sense. Thanks! nehal

Answer: 1. Results would be similar but not identical because of the stochastic nature of sequencing. I assume the same library is being sequenced at both centers or are these two separate libraries (from same sample)?
2. Nothing specific. As long as the sequencing works well data should look comparable. Technical replication for Illumina sequencing is consistent as long as libraries are appropriately handled/stored.
3. No specific effect. DNA is DNA.
4. You will get raw sequence data in fastq format of a read length equal to the number of cycles used (provided the data was not pre-trimmed). That is the standard deliverable with most sequencing centers unless you make arrangements with them to do assembly/analysis. You will get one alignment file (per sample) when the data is aligned to a standard reference. One vcf file if SNP's were called using that alignment (all chromosomes will generally be in that one file). If you ask them to do assembly then don't expect to get 46 files representing the diploid genome. It could be any number depending on how good the data was and how well the assembly worked.
5. [LiftOver][1] from UCSC.

[1]: https://genome.ucsc.edu/cgi-bin/hgLiftOver
Question: I want to dock micro RNA with a drug molecule. Please suggest some web server or online free software which can be use to study this type of docking.

Answer: RosettaDock can do it, HADDOCK can do it (free online server), many other CAPRI participants can do it.
Question: Hi.

I have usually been using hg19 reference file from ucsc.

However, ucsc has a prefix 'chr' in its files.

Of course, I am able to exclude 'chr' prefix by script programming, but I don't want to.

Instead, I need another human reference file which does not contain 'chr' prefix itself.

Do any another reference files exist not ucsc?

Answer: It's never as simple as "remove chr-prefix with a script". [hg19](http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&chromInfoPage=true) is UCSC's variant of the official [GRCh37](https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.13) assembly. Early releases of GRCh37 like [GRCh37-lite](ftp://ftp.ncbi.nlm.nih.gov:/genomes/archive/old_genbank/Eukaryotes/vertebrates_mammals/Homo_sapiens/GRCh37/special_requests/), did not use chr-prefixes, but newer releases like [GRCh37.p13](ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_genbank/Eukaryotes/vertebrates_mammals/Homo_sapiens/GRCh37.p13/seqs_for_alignment_pipelines/) adopted the chr-prefix, and use a newer mitochondrial (`MT`) sequence than hg19 does. Note also how `chrM` in hg19 is named `MT` in GRCh37. And all the unplaced contigs have very different names. So simply removing the chr-prefix in hg19 does not make it GRCh37. It makes it a wholly other chromosome naming convention, which is the last thing we need right now, since everyone has agreed to use chr-prefixes in GRCh38/hg38, except Ensembl.

**Update (Nov 4, 2016):** [Here is a UCSC Chain](https://github.com/mskcc/vcf2maf/blob/v1.6.15/data/hg19_to_GRCh37.chain) mapping UCSC's hg19 to Ensembl's GRCh37.p13 (no chr-prefix), compatible with tools like [CrossMap](http://crossmap.sourceforge.net), [Remap](https://www.ncbi.nlm.nih.gov/genome/tools/remap/docs/api), or [liftOver](https://genome.ucsc.edu/cgi-bin/hgLiftOver). Notice how all chromosomes/contigs except `chrM` only require renaming. Users of [vcf2maf](https://github.com/mskcc/vcf2maf) with hg19 VCFs as input, can pass this into the `--remap-chain` argument.
Question: Hi,

My question may seem so simple. Could you tell me what is the difference between aligning and mapping the short reads to the reference genome?

And also what is the difference between Pairwise alignment, Multiple sequence alignment and Short-Read Sequence Alignment?

Thanks

Answer: This turns out to be not so simple, so excellent question!

When we align a read, we're asking for not just where it likely came in the genome, but the exact base to base correspondence. For example, we'd like to get something like, "Read `foo` likely originated from chr1 positions 123 through 140. The first 7 bases are exact matches between `foo` and the reference, there's then a 3 base insertion, then the remaining bases match between `foo` and the reference."

When we map a read, we're just asking, "where did it come from?" We don't necessarily care about the exact alignment between the read and where it came from, though.

Until recently, "alignment" and "mapping" were pretty much synonymous. Tools like Kallisto and Salmon have changed that, since they can assign reads to genes/features/whatever without needing to look at exact alignments. Since (A) this is faster and (B) we often don't actually care about the alignment, this is a HUGE advantage in some applications.
Question: I am getting the following error in R studio when I attempt to import featureCounts count matrices into edgeR for analysis. **My question is, do I have to manually modify the output of featureCounts for use in edgeR?** Everywhere I look within documentation, it makes it seem as if I can directly load the output of featureCounts into edgeR. If the answer is as easy as switching negative values (I think featureCounts outputs '-1' in the case of no overlaps) to zero, I can handle this but it seems to me as if this is a good way to mess with statistics.

The error is:

**Error: Negative counts not allowed**

The command I am running to generate this error is:

group = c(0,1,2,3,4,5,3,4,5,3,4,5)
dge = DGEList(counts = 'file_path', group = group)

where the file path listed is the output of featureCounts run on 12 bam files. I ran featureCounts with no issue and about 60% of my reads overlapped features using the following command:

featureCounts -a 'Mus_musculus.GRCm38.95.gtf' -o features_count_all/total_file.count Sample1-1 Sample2-1 Sample3-1 Sample4-1 Sample4-2 Sample4-3 Sample5-1 Sample5-2 Sample5-3 Sample6-1 Sample6-2 Sample6-3

My matrix file looks like this. Do I need to consolidate this to be only a simple raw count number for every matrix spot rather than the comma separated fields?

ENSMUSG00000088159	1	15019040	15019159	-	120	0	0	0	0	0	0	0	0	0	0	0	0
ENSMUSG00000073737	1	15268802	15269797	-	996	2	1	2	4	1	3	2	2	2	6	2	1
ENSMUSG00000092083	1;1;1;1;1	15287254;15312363;15312452;15709485;15709485	15287484;15313030;15313030;15712548;15723750	+;+;+;+;+	15165	3	6	3	0	0	0	2	0	0	1	0	0
ENSMUSG00000102937	1	15364302	15365834	+	1533	0	0	0	0	0	0	0	0	0	0	0	0
ENSMUSG00000104149	1	15556249	15558337	+	2089	0	0	0	0	0	0	0	0	0	0	0	0
ENSMUSG00000088829	1	15685935	15686046	-	112	0	0	0	0	1	0	0	0	0	0	0	0
ENSMUSG00000077377	1	15757832	15757963	+	132	0	0	0	0	0	0	0	0	0	0	0	0
ENSMUSG00000101652	1;1	15760122;15760560	15760263;15760668	-;-	251	0	1	1	2	1	1	2	3	1	3	4

Answer: FeatureCounts adds a couple of additional gene information to the beginning of the matrix, i.e., it does not only contain the counts.
You need to remove the columns containing the gene position, strand and length (e.g. `1 15019040 15019159 - 120`) and the GeneIDs should be assigned to row names. If you read in the feature counts results first, this will become clear:

```
# read in the results (not tested, you may need to play around with read.table parameters)
fc_res <- read.table('file_path', header = T)

# assign row.names
row.names(fc_res) <- fc_res$GeneID

# exclude superfluous columns
fc_res <- fc_res[, -c(1:6)]
```
Question: Hi,

I have two matrices of normalised gene expression values from single cell data from two time points. In rows I have the same genes but cells in each matrix is not same(different time point). In figure 2 C of this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5668937/ authors say they have calculated Pearson correlation and give one Pearson correlation coefficient for each comparison as the highest correlation coefficient.

I am using many Pearson calculation like `corrplot.mixed(cor(cbind(h_160,h_140)), order="hclust", tl.col="black")` that gives me one Pearson correlation coefficient for each pair of cells (A matrix of Pearson correlation coefficient). So how I could select the highest Pearson correlation coefficient for these two matrices??

Answer: The maximum would just be:

cormatrix <- cor(cbind(h_160,h_140))
max(cormatrix, na.rm=TRUE)

Note that I have 'parallellised' the `cor()` function in R, here (see section 2 *cor function, parallellised*): https://www.biostars.org/p/273107/

In addition, it looks like they are just doing a simple network plot in Figure 2C. I also use Pearson or Spearmann correlation coefficients in my tutorial, here: https://www.biostars.org/p/285296/

Kevin
Question: Project PRJEB99111 has 147 samples. I want to download the metadata (age, sex, disease status, etc) of each sample, not fastq. The only way I can download the metadata is by downloading the xml file of each sample accession one by one - is there a way to bulk download all 147 metadata files? I can work with xml files if I have to.

You can view the metadata for a specific sample accession by clicking on the"attributes" tab. Here is an example for one sample: [https://www.ebi.ac.uk/ena/data/view/SAMEA104228123][1]

[1]: https://www.ebi.ac.uk/ena/data/view/SAMEA104228123

Answer: > Unfortunately NCBI does not contain metadata for this project.

This is not true. You can easily download a XML file containing all of the attributes of all the biosamples from NCBI. Since the procedure may also be useful in other contexts, I will describe it step by step.

First go to the page of the project (the bioproject database in NCBI speach):

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB99111

Next, get a list of all biosamples which are linked to this project. There is a section entitled "Related information" on the right site of the page. To get the list of biosamples, click on the hyperlink "Biosample".

This will open an new page which list the first 20 biosamples in the project. The URL of that page is:

https://www.ncbi.nlm.nih.gov/biosample?LinkName=bioproject_biosample_all&from_uid=400734

On the top of this page (on the right site) is a pull-down menu entitled "Send to:". Click on this menu, then select "File", then select format "Full XML (text)", and finally click on the buttom "Create File". Store the XML file on your local disk and parse it with your favorite XML tool.
Question: Hello all,

In [this paper][1] which i believe to be the original publication regarding the BigWig format, the authors write:

> Supplementary information: Supplementary byte-level details of the
> BigWig and BigBed file formats are available at Bioinformatics online.

Does anyone know where this supplementary data can be found?

Thanks and all the best :)

EDIT: None of the links for supplementary info (like the "Supplementary Table 1") seem to work in my browser. I thought this was a separate issue, but perhaps this is the real problem i'm having, and if that worked i would be able to see the Supplementary information?

[1]: https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btq351

Answer: You can find it [here][1], on the PubMed Central copy of the paper.

Oxford Press mangled most (or all) supplementary links when they moved their site to this new version.

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2922891/bin/supp_26_17_2204__index.html
Question: Hi,

I am running codeml on orthogroup tree to determine whether a certain lineage is under positive selection. The way I did is to set null model with fixed omega (`model=2`, `NSsites=0`, `fix_omega=1`) while branch-specific model by labeling branch of interest (foreground) in tree file (`model=2`, `NSsites=0`, `fix_omega=0`). After running codeml on both models, I was able to extract likelihood (lnL) and number of parameters (np) for each model, as well as Kn/Ks value for background branches and foreground branches in branch-specific model.

To found out the significant difference on Kn/Ks values between foreground and background branches, I calculate LRT by `LRT=2*(lnL1-lnL0)` and corresponding pvalues using dchisq function (df =1) in R. I noticed that in some orthogroups, LRT was reported as a negative value which leads to pvalue =0. That is definitely abnormal. How could LRT being negative happen?

Also, I am still a little bit confused about Kn/Ks value predicted for background and foreground in branch-specific model report. If I compared them and Kn/Ks is larger in foreground compared to background, may I make a comment that lineage I chose is under positive selections? The other way around, if Kn/Ks is smaller, the lineage is under negative selection.

Hope I have made my question clear. Thank you

Answer: Not completely sure I remember the details, but this may be of help. You are testing two models: one in which the omega parameter for the "foreground" branch is optimised (free parameter) and one in which omega is set at the value you choose (1, for example). The former is the alternative model and is compared to the latter, the null model, in which you fix omega at 1. The alternative model should always give you a higher or equal likelihood than the null model. The LRT test will tell you whether the difference is significant or not. In the PAML documentation there is information on how to do this test (number of parameters, one extra in this case, etc). See http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf

Regarding the interpretation of the results, no, this test does not necessarily mean whether there is positive selection or not. The LRT will tell you if the optimised omega is different from the omega you fix. It depends on the value of omega. Is it greater than 1?

I would start by freeing the omega for that branch and see if it is greater than one. If that is the case, then you want to test whether that difference compared to 1 is significant or not (calculating the likelihood with omega fixed at 1). If it is different, yes, you have evidence for positive selection. You can also test whether the omega of the foreground branch is significantly different from the general omega (model 2 vs model 0; I would start from this).

There are few cases in which the branch test will find positive selection. Positive selection is hard to detect and happens at particular sites (purifying selection at other sites may obscure the signal of positive selection). Look at the branch-site model, it has greatest power.

I hope this helps,

Federico
Question: Hello,

I'm running into "Segmentation fault (core dumped)" error when running a workflow with (or without) snakemake. I've used the same data, same tools and same parameters on several machines without a problem before.

Now the problem arises on a Ubuntu 14.04 virtual machine (4 CPU, 8Go RAM), when executing sickle. I've noticed that it depends on the parameters that I use, even when I'm specifying parameters that are supposed to be the default ones.

root@vm0083:~/mydisk/ChIP-seq_SE_GSE20870# sickle se -f fastq/GSM521934/GSM521934.fastq -t sanger -o sickle.fastq -q 20 -l 20

FastQ records kept: 2449731
FastQ records discarded: 1137314

root@vm0083:~/mydisk/ChIP-seq_SE_GSE20870# sickle se -f fastq/GSM521934/GSM521934.fastq -t sanger -o sickle.fastq -q 20 --length-threshold 20
Segmentation fault (core dumped)

root@vm0083:~/mydisk/ChIP-seq_SE_GSE20870# sickle se -f fastq/GSM521934/GSM521934.fastq -t sanger -o sickle.fastq --qual-threshold 20 -l 20
Segmentation fault (core dumped)

Parameters names and default values:

root@vm0083:~/mydisk/ChIP-seq_SE_GSE20870# sickle se -h
sickle: invalid option -- 'h'

Usage: sickle se [options] -f -t -o

Options:
-f, --fastq-file, Input fastq file (required)
-t, --qual-type, Type of quality values (solexa (CASAVA < 1.3), illumina (CASAVA 1.3 to 1.7), sanger (which is CASAVA >= 1.8)) (required)
-o, --output-file, Output trimmed fastq file (required)
-q, --qual-threshold, Threshold for trimming based on average quality in a window. Default 20.
-l, --length-threshold, Threshold to keep a read based on length after trimming. Default 20.
-x, --no-fiveprime, Don't do five prime trimming.
-n, --trunc-n, Truncate sequences at position of first N.
-g, --gzip-output, Output gzipped files.
--quiet, Don't print out any trimming information
--help, display this help and exit
--version, output version information and exit

I would like to know if this issue is likely related to sickle itself, or the VM I'm using, or something else...

Thank you

Answer: In case someone else is running into this issue, using shortcuts for params names avoids the problem: `-q` instead of `--qual-threshold` and `-l` instead of `--length-threshold`, specifically.

So in my case the problem was not related to my VM, but to the fact that my VM was using another version of my code which included full params names.
Question: I have some paired end short read sequencing data. I'd like to know the fragment alignment overlap between the paired reads, and prune that part out for one of the reads to avoid processing the same region twice. For my specific downstream purposes, it has to be done in C.

My current idea is to use the `MC` tag, which contains the CIGAR string of the mate read, and pass it to `bam_cigar2rlen()`. Together with `core.mpos`, I should be able to determine the end position of the mate read.

However, `bam_cigar2rlen()` takes `int n_cigar` and `const uint32_t *cigar` for the arguments, and it is unclear to me how to get these parameters from the `MC` tag. I am also not sure if parsing the CIGAR string like this is the most efficient way.

I'm a novice to C and I'd highly appreciate any help. Thanks!

Answer: Unless your data files have another (locally-defined) tagged field giving the mate's end position directly, computing it by parsing the MC tagged field and adding the length returned by `bam_cigar2rlen()` to MPOS is indeed the best way to compute this.

As you've seen, `bam_cigar2rlen()` uses HTSlib's already-parsed internal CIGAR representation. Since v1.12, HTSlib has `sam_parse_cigar()` for parsing a C string (as returned by `bam_aux2Z()` on your MC field) into the representation you need.
Question: Hi,
I have a set of ~150 gene and I want to obtain a list of phenotypes (diseases or conditions) they are associated with. I was wondering what is the way to do this in OMIM in a batch manner. Eg.,

Gene Condition

BRCA1 Breast Cancer

Thanks for your help!

Answer: If you are willing to use R and biomaRt:

gene_mart <- biomaRt::useEnsembl(biomart = "ENSEMBL_MART_ENSEMBL",
host = "www.ensembl.org",
dataset = "hsapiens_gene_ensembl")

biomaRt::getBM(mart = gene_mart,
attributes = c("hgnc_symbol", "mim_morbid_accession", "mim_morbid_description"),
filters = "hgnc_symbol",
values = c("BRCA1"),
uniqueRows = TRUE)

This will give the following output:

hgnc_symbol mim_morbid_accession mim_morbid_description
BRCA1 114480 BREAST CANCER;;BREAST CANCER, FAMILIALBREAST CANCER, FAMILIAL MALE, INCLUDED
BRCA1 604370 BREAST-OVARIAN CANCER, FAMILIAL, SUSCEPTIBILITY TO, 1; BROVCA1BREAST CANCER, FAMILIAL, SUSCEPTIBILITY TO, 1, INCLUDED;;OVARIAN CANCER, FAMILIAL, SUSCEPTIBILITY TO, 1, INCLUDED
BRCA1 614320 PANCREATIC CANCER, SUSCEPTIBILITY TO, 4; PNCA4
BRCA1 617883 FANCONI ANEMIA, COMPLEMENTATION GROUP S; FANCS

If you want more than one gene, you just have to write all your hgnc symbols in the *values* parameter, example with 2 genes:

biomaRt::getBM(mart = gene_mart,
attributes = c("hgnc_symbol", "mim_morbid_accession", "mim_morbid_description"),
filters = "hgnc_symbol",
values = c("BRCA1", "FOXP3"),
uniqueRows = TRUE)

This script will get the omim IDs and descriptions filtered by hgnc symbols. If you prefer to use ensembl gene IDs you can replace *filters = "hgnc_symbol"* by *filters = "ensembl_gene_id"*. For more information you can refer to the biomaRt user guide at https://www.bioconductor.org/packages/devel/bioc/vignettes/biomaRt/inst/doc/biomaRt.html
Question: How to merge unique/non overlapping genes between 2 gene model GFF3 files?

I have 2 gene models:

1. Recent supposedly higher quality gene model
2. Older supposedly lower quality gene model

Some genes though that are described in literature have been removed in the more recent gene model. While they are in the older gene model.
I would like to add these and other "missing" genes back to the more recent gene model, from the older gene model.

But if there is an overlap between genes from the different models, I would like to keep the genes as described in the more recent gene model.

Is there a way to do this?

Answer: `agat_sp_complement_annotations.pl` from [AGAT][1] does the job also. It allows different coding gene loci to overlap at their UTR levels ( because it is often less well annotated ). It allows also to get e.g a non-coding element ( ncRNA, pseudogene, repeat) to overlap a coding gene. So it a bit less strict than bedtools.

[1]: https://github.com/NBISweden/AGAT
Question: I want to know if there is a way to keep sites with 0% missing a VCF file.

On vcftools I can do:

`vcftools --vcf input_file.vcf --max-missing 1.0 --out output_noMissing`

Is there a way to do this with BCFTOOLS as it is much faster?

Answer: not tested:

bcftools view -g ^miss input.vcf
Question: I got a report from our vendor for the RNA samples for the RNA seq analysis we sent.
one of the samples from a group has a RINe score of 5.3. Should this sample be resent before we sequence them?

Answer: The consensus in the scientific community is a minimum RIN of 8. Samples with RIN < 7 don't work well in RNA-seq or other NGS. I would redo that sample.

Best,

Rafael
Question: Lets say i have RNA sequence data from three groups: \
1 = Control_untreated (healthy cells), \
2 = Case_untreated (diseased cells), \
3 = Case_treated (diseased cells with siRNA knockdown against a gene suspected of causing disease by altering the expression profile)

Question 1:
Would i be able to use a PCA plot to summarise the RNA expression profiles and see if the treatment drags the case group towards the control group?

Question 2 (if Q1=yes):
I know oxford nanopore systems are not ideal for RNA-seq yet due to lacking dept. But since im interested in comparing just a summery of the expression profiles would it be possible to use it for generating the data i need?

Question 3 (if Q1=yes):
Is it possible to do something similar for an DEU or DTU analysis?

Answer: 1. Yes you can. You can also try TSNE or ICA
2. If you just care about expression, nanopore is not a good choice. Traditional RNA-seq or microarray is better. Nanopore is designed for long-read sequencing and has advantages in calling variants such as splicing isoform, gene fusion. Quantifying the expression level is not the aim of nanopore.
3. Yes