cannin/biostars_qa huggingface dataset:
Question: I want to find the coordinates of all occurrences of the sequence recognized by a restriction enzyme. I know that using EMBOSS I may do this, but this task seems perfectly fitted for short-read sequence alignment software. However, I didn't find any reference for this.

I used bwa for the task and quickly obtained some results. However, to be on the safe side I will like to ask is someone has done something similar or has some advice, perhaps I am stretching the use of bwa.

I tried the following:

echo -e ">DpnII\bGATC" > DnpII.fa
bwa aln -N -n 0 -o 0 -e 0 -l 4 -k 0 dm3 DpnII.fa > DpnII.sai
bwa samse -n 100000000 dm3 DpnII.sai DpnII.fa > DpnII.sam

The results seem to match the right sites, also the number of sites (489570 for Drosophila) that I obtain are close to what I expect.

Thanks.

Answer: Here is an update for people landing on this question: The bwa method that I proposed does not work properly. I don't remember exactly why it was wrong but I stop using it. I think it was returning regions with NNNs.

My solution was to use Biopython to search for a pattern. This is quite fast and accurate. Here is the code that finds and sorts the results:

https://github.com/maxplanck-ie/HiCExplorer/blob/master/hicexplorer/findRestSite.py

[update]
Based on the comments by [Michael Dondrup][1] I fixed an error in the code.

[1]: https://www.biostars.org/u/55/
Question: Hello all,

I am looking at the Level 3 CNV files on TCGA - the ones generated using SNP microarrays. I have a few questions:

1. How is 'segment mean' calculated and what is the exact biological interpretation?
2. For each patient I have two files called e.g. `....hg19.seg.txt` and `...nocnv_hg19.seg.txt`. What does each file contain, and which should I be using?

Thanks for any help,

Stephanie

Answer: [Here][1] is the best documentation (that I know of) for TCGA SNP-array based CNV data. Regd. your two questions:

1. CBS segmentation algorithm identifies regions in the genome that, in spite of noise, probably have a uniform underlying copy number. The "segment mean" of each region is reported in the level 3 file, and can be used as the estimated CN-ratio for the segment.
2. "nocnv" just means that germline CN variations are removed. In TCGA, they ran SNP arrays on normal tissue too.

[1]: http://software.broadinstitute.org/cancer/software/genepattern/affymetrix-snp6-copy-number-inference-pipeline
Question: Hi,

I wanted to download 10 abstracts for 10 specific mouse genes (for a total of 100 abstracts). I was thinking of searching Pubmed for it, so I put in my first gene (FOS) in the search bar. But there was no option (that I could find) for mentioning the species I want to restrict my search to. So the results were either for 'Humans' or 'All organisms'.

Since this is a pretty basic feature, I was hoping that Pubmed would allow it inherently. But it wasn't even mentioned in the 'Advanced' options. So, does anyone know if and how this can be done?

Answer: This thread piqued my interest since I wanted to write a script to do this for a web app. I'm not sure why the above examples don't work but here's how I got it going.
First you need the gene ID for Fos. You can do that just by searching:

https://www.ncbi.nlm.nih.gov/gene/14281

The Gene ID is the number at the end of the URL above.
You can put this directly into pubmed as the `&from_uid=` in this url:

https://www.ncbi.nlm.nih.gov/pubmed?LinkName=gene_pubmed&from_uid=14281

You can also use elink as suggested by Pierre. To do that just tack it onto this url on the `&id=`:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=gene&dbto=pubmed&cmd=neighbor&retmode=xml&id=14281

Change 'gene' to 'protein' if you use protein accession numbers. This will return an XML like so:

gene

14281

pubmed
gene_pubmed

26762887

26143639

etc...

The pub med hits are in ``. You can extract these with a parsing tool or do it by hand for a few genes. Then put them into pubmed as a unique id:

https://www.ncbi.nlm.nih.gov/pubmed/?term=26762887

To put this in a program you could loop the gene names and paste them into either of the links above. For the eutils link you need to capture the url and parse it, extract the ids and paste them into a url.
Question: I've got a bam file of paired-end illumina sequences mapped to a reference sequence and I was wondering if when I call the samtools view command if the results represent only reads that completely span over the region of interest or if it also included reads that intersect it by one or two base pairs.

If not are there any existing tools out there to *only* read thats completely span over the region of interest?

with this visual example of what I mean I would want to only include sequences like case_2 and exclude case_1:

[===================Region-of-Interest=====================]

[---------------------------------------------------------------------------------------------------------] case_1
[-------------------------------------------------------------------------------------------------------------------] case_2

Thanks!

Answer: [bedtools intersect][1] with -f 1.0 or -F 1.0 should do what you want. Something like:

bedtools intersect -F 1.0 -a intervarls.bed -b file.bam

I am not an expert on bedtools, though, so you may have to tweak / correct this command.

[1]: http://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
Question: Hello,

I am looking at the human genome available through Bioconductor packages: NCBI GRCh38 and UCSC.hg19. And I do not get all the different sequence names I see. Could you help please?

In the UCSC.hg19, I do have chromosomes 1 to 22 + X and Y. But I also have chrM, chr1_gl000191_random, chr4_ctg9_hap1, chrUn_gl000212....and so on.

In the NCBI GRCh38, I can see sequences called MT, HSCHR1_CTG2_UNLOCALIZED, HSCHR3UN_CTG2, HSCHR2_RANDOM_CTG1....

What are those _random, _unlocalized, chrUn, _ctg9_hap1...please??

Should all those sequences be used when trying to align NGS reads to the genome for instance? or only a subset?

many thanks

Aurelie

Answer: chrM == MT == mitochondrial DNA. You should probably use this.

The `chrUn_*` sequences are unplaced contigs. So they may belong in the genome, but we don't know where. I personally use these, but I know that's not universal.

`chr*_unlocalized` and `chr??_*_random` are contigs that are known to belong to a specific chromosome (I've never looked into how that was determined) but haven't yet been integrated in. You'll want to use these.

The various `*hap*` chromosomes are alternate haplotypes. There aren't a lot of good ways to deal with these during alignment yet, so a lot of people don't use these. An upcoming version of BWA is supposed to handle these in a good way, see if you plan to use BWA keep an eye out for that update and then definitely use these.
Question: Hello everyone,
I would like to know if there are tools that feeded with ChIP-seq data for various histone marks, identify the enhancers by the "right" methylation/acetylation combinations.
Until now, I found CSI-ANN, but I was wondering if there are other tools.
Thanks

Answer: The search term you're looking for is "chromatin segmentation". For example, omicstools has [a list of chromatin segmentation tools][1]

Of these, [ChromHMM][2] is probably the most widely used in publications at the moment.

Also see other answers on biostars on chromatin segmentation and enhancers:

- https://www.biostars.org/p/185468
- https://www.biostars.org/p/17912
- https://www.biostars.org/p/175982
- https://www.biostars.org/p/142111

[1]: https://omictools.com/chromatin-segmentation-category
[2]: http://compbio.mit.edu/ChromHMM/
Question: Hi everyone,

I would be very grateful if you could help me.

I want to download the sequences of all the exons for each human gene. I went to ensembl biomart and tried to do it for BRCA2 first (this was a random choice). First I selected the following attributes: Unspliced(gene), Exon start, Exon end, strand, Gene start. I thought that this information will be enough to 'cut out' all the exons. The first thing I noticed is that the number of exons is almost twice larger then the number I see on the Wiki.
Then I tried to download directly Exon sequences, but the number of sequences was again larger the the number of exons should be and moreover several exon ids correspond to the same sequence.
I also tried to download coding sequence and cDNA, but the length of both sequences is not consistent with 'official' BRCA2 length!

These all drives me crazy and I have absolutely no idea on what to do. All I want is to get for each gene the sequences of its exons. Help me please!

Thanks!

Answer: Download refGene.txt.gz from UCSC :

wget http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/refGene.txt.gz

Then keep only uniq Gene name ( column 13 ) and extract coordinate ( chrom - exonstart - exonend) to a bed file. Column 10 and 11 contains exonsStart and exonEnd position separated by a comma.

zcat refGene.txt.gz|sort -u -k13,13|cut -f3,10,11|awk 'BEGIN{OFS="\t"}{split($2,start,",");split($3,end,","); for(i=1;i exons.bed

Finally, you can get exons sequence using bedtools getFasta :

bedtools getfasta -fi hg19.fa -bed exons.bed -fo exon.fa

paid attention to work only with chromsome name in range 1-22,X,Y.
Question: Hello All,

I am confused about RA-seq normalization methods and when is appropriate to use them. I appreciate if you could share your thought with me.

My understanding from TMM and TPM is that TMM is appropriate for between sample/condition comparison as it counts for RNA composition in addition to library size (e.g. it is used by edgeR for DE analysis). While a method like TPM is better for within sample comparison. But many online tools use TPM to illustrate a gene expression levels across different tissue types (like GTEx data portal). Isn't this wrong?

The reason that I got into this is that I was performing DE analysis using edgeR the other day to compare samples from a tumor to some normal tissues. Although edgeR reports a downregulation for some genes, an illustration of those gene's expression using TPM values shows kinda an upregulation effect. When I extracted and plotted the pseudo counts from edgeR (TMM normalized counts), I clearly can see the downregulation - but TPM doesn't agree! So I am confused now. Is it wrong to use TPM normalized counts for plotting gene expression?

I appreciate your time and thoughts on this.
Thanks,

Answer: Online tools use TPMs for illustration because they can calculate them once and be done. This is convenient when people want to add samples over time and change groups and samples being visually compared. The results of that will not be as robust to outliers as normalized counts (produced with TMM or another method), but they're usually good enough for visualizations.

What you observed is due to the non-robustness of TPM that I mentioned earlier. TPM is one of those things that has its use, but if you're in a scenario where you can use properly normalized counts then that's usually preferable. As an aside, you could convert your normalized counts (or pseudo-counts) to TPMs and then you'd see the down-regulation.
Question: Dear community,

I ran phyml on a gene family to build a tree. Looking at the results, I'm a bit worried about the log-likelihood value: it's -754, which means the likelihood is almost zero! Does this mean that the program has little confidence in the estimated parameters or the tree topology? I was wondering if I'm understanding this incorrectly.

Thank you so much!!

```

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
--- PhyML 3.3.20190909 ---
http://www.atgc-montpellier.fr/phyml
Copyright CNRS - Universite Montpellier
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

. Sequence filename:	exon3_wb_aligned_phy
. Data set:	#1
. Initial tree:	BioNJ
. Model of nucleotides substitution:	GTR
. Number of taxa:	52
. Log-likelihood:	-754.13849
. Unconstrained log-likelihood:	-348.41766
. Composite log-likelihood:	-6438.01266
. Parsimony:	119
. Tree size:	1.35199
. Discrete gamma model:	Yes
- Number of classes:	4
- Gamma shape parameter:	1.901
- Relative rate in class 1:	0.28116 [freq=0.250000]
- Relative rate in class 2:	0.64406 [freq=0.250000]
- Relative rate in class 3:	1.06730 [freq=0.250000]
- Relative rate in class 4:	2.00748 [freq=0.250000]
. Nucleotides frequencies:
- f(A)= 0.37232
- f(C)= 0.24092
- f(G)= 0.17327
- f(T)= 0.21350
. GTR relative rate parameters :
A <-> C 0.82212
A <-> G 1.82689
A <-> T 0.53724
C <-> G 0.17829
C <-> T 2.00016
G <-> T 1.00000
. Instantaneous rate matrix :
[A---------C---------G---------T------]
-0.82453 0.25951 0.41474 0.15028
0.40104 -1.00102 0.04048 0.55950
0.89119 0.05628 -1.22720 0.27973
0.26208 0.63136 0.22702 -1.12046

. Run ID:	none
. Random seed:	1625516914
. Subtree patterns aliasing:	no
. Version:	3.3.20190909
. Time used:	0h0m4s (4 seconds)

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
Suggested citations:
S. Guindon, JF. Dufayard, V. Lefort, M. Anisimova, W. Hordijk, O. Gascuel
"New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0."
Systematic Biology. 2010. 59(3):307-321.

S. Guindon & O. Gascuel
"A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood"
Systematic Biology. 2003. 52(5):696-704.
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooooooo

```

Answer: The principle of maximum likelihood is to choose the tree which makes the data most probable. As tree probabilities are usually tiny, especially for large datasets, we express them as ln(P), which is the log likelihood (LL). LL is a negative number (log function is negative in the 0-1 range), and the best it can be is 0 (when the probability is 1).

Seems like you have a smallish tree, and the LL you obtained is appropriate. As of this writing I am monitoring an ongoing large tree that has `LL=-978001.783`, so you are golden. LL values are comparable for the same dataset (alignment), but not between different datasets. Higher LL is better.
Question: Hi community!,

I'm annotating variants with the VEP software and I'm finding some unexpected transcript data of the type:

- NM_014938.4_dupl16
- NM_001170637.2_dupl3

```
1 206516261 . C T 47 PASS CSQ=T|non_coding_transcript_exon_variant|MODIFIER|SRGAP2|23380|Transcript|NM_001170637.2_dupl3|mRNA|1/20||NM_001170637.2_dupl3.1:n.65C>T||65|||||||1||SNV|EntrezGene|||||||||||C|C||||||||||||||||||||||||||||||||||||||||||||||||||||||||||1:206516261-206516261|0.4996565||||,
T|missense_variant|MODERATE|SRGAP2|23380|Transcript|NM_001170637.3|protein_coding|1/20||NM_001170637.3:c.65C>T|NP_001164108.1:p.Arg289Trp|864|865|289|R/W|Cgg/Tgg|||1||SNV|EntrezGene||||||NP_001164108.1|||||C|C|OK|||||||||||||||||||||||||||||||||||0.63580||||T|T||||||||||2|||||||1:206516261-206516261|0.4996565||||,
T|missense_variant|MODERATE|SRGAP2|23380|Transcript|NM_001300952.1|protein_coding|1/18||NM_001300952.1:c.65C>T|NP_001287881.1:p.Arg289Trp|864|865|289|R/W|Cgg/Tgg|||1||SNV|EntrezGene||||||NP_001287881.1|||||C|C|OK|||||||||||||||||||||||||||||||||||0.63580||||T|T||||||||||2|||||||1:206516261-206516261|0.4996565||||,
T|non_coding_transcript_exon_variant|MODIFIER|SRGAP2|23380|Transcript|NM_015326.3_dupl3|mRNA|1/20||NM_015326.3_dupl3.1:n.65C>T||65|||||||1||SNV|EntrezGene||YES|||||||||C|C||||||||||||||||||||||||||||||||||||||||||||||||||||||||||1:206516261-206516261|0.4996565||||,
T|missense_variant|MODERATE|SRGAP2|23380|Transcript|NM_015326.4|protein_coding|1/20||NM_015326.4:c.65C>T|NP_056141.2:p.Arg289Trp|864|865|289|R/W|Cgg/Tgg|||1||SNV|EntrezGene||YES||||NP_056141.2|||||C|C|OK|||||||||||||||||||||||||||||||||||0.63580||||T|T||||||||||2|||||||1:206516261-206516261|0.4996565|||| GT:DP:VD:AD:AF:RD:ALD 0/1:9:3:6,3:0.3333:6,0:3,0
```

Searching on the VEP webpage or in the internet I can't find any reference to this kind of "dupl" suffix. Has anyone faced this? I don't know if they are alternatives of the transcript or explain why they are not transcripts on its own.

Thanks in advance!

Cristian.

Edit: Added example of variant with the vep annotation of dup (NM_015326.3_dupl3)

Edit2: Using VEP ensembl version 91.1 with cache v91

Answer: We are investigating these. It looks like some RefSeq transcripts (eg NM_001170637.3) have been duplicated in Ensembl's other_features database with a lower version number and this dupl suffix (eg NM_001170637.2_dupl3). This has been propagated across to the VEP cache, which is why you're seeing them. We don't currently know why, but we believe that you can just ignore them from your analyses for now.
Question: One of our project used to query OMIM data as XML through NCBI's efetch utility, as described here for example:

What is the best way to interact programmatically with OMIM?

However, it seems the service has stopped functioning a few months ago. It now simply returns the following error:

Database: omim - is not supported

I can find no mention of an update to the API on NCBI's website or anywhere else.
At the same time, the pages accessible directly on OMIM's website offer no link to structured data (XML or otherwise) and the downloadable file, while using some specific format to delimit fields, is still far from the flexibility of the former XML files (for example, it is impossible to retrieve metadata for each reference).

Is there currently any way to regain access to OMIM data in a structured, parsable format (XML...)?

Answer: FYI: this was just announced this morning via twitter:

https://twitter.com/#!/OmimOrg/status/196939511220015104

@OmimOrg

OMIM API is now open, see http://omim.org/api and
http://omim.org/help/api #OmimOrg
Question: Hi guys,

I have an R programming question. I want to compare the genotype (`.GT` columns) with the given alleles (`.allele` columns) in dataframe `df1` and see if they are concordant or not. The rule is that if there is only one allele, it should be 0/0. If there are two alleles (for example GA), the genotype should be 0/1 and that is why I have mismatch in the concordance column for GA. So, A, T, G or C individually is 0/0 and in pair combination is 0/1. Based on this information, I want to add new concordance column next to every two columns that are compared with and have match or mismatch result in concordance column. This concordance column should be cbind next to every compared pair of column and get `Result`. Could you please help me get this done. Thank you.

`df1`

```
1.allele 1.GT 2.allele 2.GT
A 0/0 A 0/0
GA 0/0 CT 0/1
C 0/0 G 0/0
```

`Result`

```
1.allele 1.GT 1.Concordance 2.allele 2.GT 2.Concordance
A 0/0 match A 0/0 match
GA 0/0 mismatch CT 0/1 match
C 0/0 match G 0/0 match
```

Answer: Data:

```r
df <- data.frame(allele1=c("A","AT","C"),GT1=c("0/0","0/0","0/0"),
allele2=c("AT","G","CG"),GT2=c("0/0","0/0","0/0"))
```

First I create a vector of correct translation. But I am assuming that there is always a 0/1, never a 1/0:

```r
library(stringr)
translation <- function(x) ifelse(str_length(df[,x])>1,"0/1","0/0")
correct <- sapply(c(1,3),translation)
```

Finally, mis/match:

```r
match.fun <- function(x) ifelse(df[,x] == correct[,(x/2)] , "match" , "mismatch" )
comparison <- sapply(c(2,4),match.fun)
cbind(df[,1:2],"Concordance1"=comparison[,1],df[,3:4],"Concordance2"=comparison[,2])
```
Question: Dear all

I followed some links here in biostar to get the differential expressions of my RNAseq data for tumor vs control.
Then I get the pathways, I did somatic mutations using GATK pipeline to get some somatic mutations.

I found some differentially expressed genes and found common somatic mutations in them, could be interesting.
Then I analyzed the top pathways to see if they are related to cancers, nothing interesting is found.

I am still trying to connect pieces. Any suggestion how can I conclude my results? What else we can do?

Thank you

Answer: Just some ideas off the top of my head:

1. **Mutation-to-expression modelling:** For each mutation, test it's association to the expression of
differentially expressed genes (DEGs) in the mutation's 'vicinity'.
This can be as easy as building a linear regression model with
expression as the *y* (dependent) variable and mutation
`present`/`absent` as *x* (predictor). From this, you could derive
R-squared values and cross validated 'shrunk' R-squared values,
along with p-values. y variable would be continuous; x variable
would be categorical with mutation absent as reference/base level.
2. **Transcription factor binding sites:** Check for new TFBS (transcription factor binding sites) that may be
introduced as a result of each mutation. Look at databases like
JASPAR to do this - there are also other threads on biostars. There
are undoubtedly some mutations in your data that are going to
modulate expression of nearby genes. For an idea of mechanism, see
the wonderful study by Manour: An
Oncogenic Super-Enhancer Formed Through Somatic Mutation of a
Noncoding Intergenic Element
3. **Histone binding regions**: Check for overlapping histone methylation (e.g. H3K27me3) and
acetylation (e.g. H3K27ac) binding regions - this data is available
from the UCSC, as far as I know. A mutation in such regions could
modify chromatin structure and alter expression.
4. **Transcription start sites:** Overlapping transcription start sites (TSS) - again,
available from UCSC I believe
5. ***In silico* prediction:** Use one of those functional / pathogenicity prediction tools. There
have been many tools released in recent years, including ones
tailored for cancer and somatic mutations. Take a quick look here:
https://www.biostars.org/p/286364/#286483

Noe that, technically, you could introduce all of the data from points 2-5 into the model mentioned in point 1. This would then be a robust way to assess the role of each mutation in relation to gene expression.

Finally, thinking just about the RNA-seq data, you could deconvolute it in order to identify immune cell-types that may be present in the tumour. This would give you an indication of the amount of immune cell infiltration, which is likely to differ across your tumors.

There are yet more ideas that I have not mentioned.

Kevin
Question: scRNA-seq novice here: We have four 10X scRNA-seq samples (wildtype and knockout condition) as n=2 each.
Each pair (so one WT and one KO) was produced on the same day respectively, same FACS sorting machine, same lab, same technician etc, so avoiding batch effects as much as we could.

For comparative analysis between the conditions I went through the `scran` / [OSCA][1] workflow and now aim to integrate the datasets. Essentially the choice is now to either merge the datasets without explicit batch correction via fastMNN (and only do per-sample depth correction via `multiBatchNorm` to ensure equal depth across the already normalized samples) or to apply fastMNN. I tested and visualized both approaches for every replicate independently, see below, and see quite different results.

Both replicates (if no fastMNN is applied) show a reproducible separation by condition (which we expect), so probably the influence of condition is greater than any batch effect. When applying fastMNN the two conditions lose this separation.

**Therefore my question**: Are there situations where batch correction masks interesting biological features. Given that we see reproducible separation by condition, could it be more meaningful to not apply fastMNN? If I combine the datasets and only correct cor batch = day (so rep1 is one batch and rep2 is one batch) I manage to preserve the separation by condition. The tSNEs then pretty much look like the left panel in the plot below.

Comments and your experiences with this are appreciated.

![enter image description here][2]

[1]: https://osca.bioconductor.org/
[2]: https://i.ibb.co/M1kf9xM/Rplot.png

Answer: I would apply the mutual nearest neighbours correction for the exact reasons laid out by the developer here: https://osca.bioconductor.org/multi-sample-comparisons.html#sacrificing-differences