cannin/biostars_qa huggingface dataset:
Question: Hi,

I was wondering which NCBI reference genome assembly to use for mouse GRCm38, if I don't want to use the UCSC mm10. I thought the FTP-site of the Sanger mouse genomes project might be a good place to check: ftp://ftp-mouse.sanger.ac.uk/ref/ Does anyone know what the 68 refers to in the file name - [GRCm38_68.fa][1]?

Many thanks,
Lorna

[1]: ftp://ftp-mouse.sanger.ac.uk/ref/GRCm38_68.fa

Answer: This is one, but from ensembl --> ftp://ftp.ensembl.org/pub/release-68/fasta/mus_musculus/dna/. Inappropriate?

This Ensembl version (release 68) corresponds to GRCm38 (mm10) (https://www.biostars.org/p/52254/)

Some sorces to download mm10 genome version:

http://hgdownload-test.cse.ucsc.edu/goldenPath/mm10/bigZips/

http://www.ncbi.nlm.nih.gov/genome/52
Question: I am new to analysing ATAC-seq data.

As mentioned on https://www.biostars.org/p/209592/, there seem to be two ways to use MACS to analyse ATAC-seq data.

1. Utilising the --shift -100 --extsize 200 command
This, I believe, is to find where the cutting sites are.

2. Utilising the --shift 37 --extsize 73 command
This is to find the nucleosomes since the DNA is wrapped around nucleosomes is circa 147bp.

I specifically want to document regions of open chromatin, particularly enhancer regions.

I have read in the literature that open chromatin regions will give reads <100bp. If this is correct, do I have to filter my bam file to look at only this length.

Surely if I map all of the my data, not filtering for read size I will end up with genome-wide ATAC-coverage? Looking at my bam file following mapping to the genome, my read sizes range from 50bp-450bp (with a peak at around 70bp).

Answer: You can use MACS, but it is not the best tool for identifying open chromatin regions. I would recommend using Fseq: http://fureylab.web.unc.edu/software/fseq/ It was designed to analyze DNase -seq data, the predecessor of ATAC-seq. When analyzing data from either experiment, you do not really want to look for overlapping reads the way you do for Chip-seq. You want to look for regions with a high density of cut sights or insertion sights, and Fseq is designed to highlight those.
Question: Hi,

By default IGV indicates insertions with respect to the reference with a purple I, and I would like to hide it or change the color, I have read the manual and it is making me crazy, I can´t find the way to do it, somebody knows if it is possible?

Answer: Preferences -> Alignments -> Hide indels < 'x' bases (?)
Question: Hi,

I am trying to compare few models for evoultionary distance calculation. I want to compare just those basic ones - Jukes-Cantor, K2P, Tamura, Tamura-Nei, HKY(Hasegawa-Kishino-Yano), GTR (general time reversible = Tavare) and Falsenstein81. I have already found equations for evolutionary distance calculations for first 4 models mentioned above. I am still looking for HKY,GTR,F81. Can you help me? I mean I know how to calculate substitution rates (for transitions and transversions) but please help me how to calculate evolutionary distance then?

Thanks Sam :)

Answer: I would first look at the original papers for the different models, you will find the equations for sure. I know that the APE/phangorm package uses those but I think they merely call [phyml][1]. I looked briefly and you can request the source code and that will give you some code to work with.

[1]: http://www.atgc-montpellier.fr/phyml
Question: Hi guys,

I've to create a variable a in which elements are column names of a database.

For example, I have the following column names of a df:

colnames(df)
"PS_01", "PS_01_mod2", "PS_02", "PS_02_mod2"

I want to create a vector a in which elements are the column names of df ending with mod2, so I want to create this situation

a <- "PS_01_mod2", "PS_02_mod2"

How can I do this in a simple way?

Answer: Hi, a simple `grep()` command will do this for you, and we can add a regular expression ('regex') to ensure positional specificity:

vector <- c('PS_01', 'PS_01_mod2', 'PS_02', 'PS_02_mod2', 'mod2_mod1')
idx <- grep('mod2$', vector)
vector[idx]
[1] "PS_01_mod2" "PS_02_mod2"

The dollar, `$`, means that we only want '*mod2*' appearing at the end of a line.

Note the difference here, without the dollar:

idx <- grep('mod2', vector)
vector[idx]
[1] "PS_01_mod2" "PS_02_mod2" "mod2_mod1"

Kevin
Question: Hi!

I have an alignment (.bam) of reads to mm9 genome. I sorted it with `samtools sort`, so that later I can use `-sorted` key with bedtools. I also created a .bed-file with regions of interest, in which I want to count number of reads, that mapped to them. I tried this: converted .bam to .bed with `bedtools bamtobed`, and then intersected them counting number of hits (`bedtools intersect -a regions_of_interest.bed -b alignment_sorted.bed -c -sorted > Neg2H_counts.bedgraph`). The problem is, it looks fine for all chromosomes with numbers from 0 to 9 (and X), but all counts for all regions of interest of chromosomes with higher number (chr10, chr11, etc) are 0. There is no biological reason for that, in fact the highest signal should be on chr11. What could be wrong here? I am fairly new to all these tools.

**UPDATE**

I tried to do the same intersection with bedmap and the result is identical... So there probably is something wrong with my files - what could it be?

I also tried sorting the alignment-derived bed-file in the same way, as I did with the files with regions of interest and it doesn't help.

Answer: By default, the `-sorted` option assumes that both of your input files have chromosomes sorted lexicographically (chr1, chr10, chr11, etc.). I suspect one of your files is this way and one is not. Alternatively, if you want to use a different order, you can use the `-g` option. See the docs for details: http://bedtools.readthedocs.org/en/latest/content/tools/intersect.html#g-define-an-alternate-chromosome-sort-order-via-a-genome-file
Question: This is a follow-up to an earlier question of mine.

I tried the following command (copied verbatim from [UCSC's mysql access docs][1]):

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A

...but it fails with the error message:

ERROR 1045 (28000): Access denied for user 'genome'@'my.host.redacted' (using password: YES)

What else must I do to get through?

Thanks!

[1]: https://genome.ucsc.edu/goldenPath/help/mysql.html

Answer: OK, problem solved.

**Diagnosis:** PEBKAC [1]

**(one) Solution:** change the command line to

mysql **--no-defaults** --user=genome --host=genome-mysql.cse.ucsc.edu -A

---

[1] More specifically: I had a long-forgotten `~/.my.cnf` file lingering around, which was required by a pipeline that I used briefly back in 2006; even adding `--verbose` to the command line did not reveal that an incorrect connection parameter, namely a password, was being read from this file; I knew, from the remark `"(using password: YES)"` in the error message, that a password was somehow being used, but I had no way of determining the source of this password, nor what it was. How I finally figured out the problem is a tl;dr-length saga I won't bore you with.
Question: Hello, I always believed expression matrix on GEO is normalized. However, I get huge big log2FC from [GSE85957](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85957) today.

> head(expr_3)
# A tibble: 6 x 8
SYMBOL logFC AveExpr t P.Value adj.P.Val B ENTREZID

1 Spp1 3419. 2421. 10.2 0.00000368 0.00128 0.221 25353
2 Gstp1 2125. 2338. 10.3 0.00000328 0.00128 0.243 24426
3 Cyp2e1 2047. 2833. 4.32 0.00204 0.0235 -1.89 25086

Here is how my expression data extracted

gse_path <- "/datapool/pengguoyu/Microarray/20190711_geo/rawdata/GSE85957_series_matrix.txt.gz"
gse <- getGEO(filename=gse_path, AnnotGPL=TRUE)
expr <- exprs(gse)
So I go back to check expression matrix

PROBEID GSM2288460 GSM2288461 GSM2288462 GSM2288463 GSM2288464
1367452_at 1165.0328 1011.4838 1193.8429 1143.6874 1162.2721
1367453_at 512.07166 519.57355 502.8087 433.26254 480.2318
1367454_at 647.18243 619.50635 673.89526 644.89575 685.5907
1367455_at 1226.1555 1299.9249 1318.0239 1363.5055 1308.6063
1367456_at 1530.6841 1611.0748 1768.4469 1761.0474 1751.5911
1367457_at 426.08826 282.9359 433.74475 421.27148 445.81595

This seems to be data without normalized. How can I know any one expression matrix from GEO is normalized or not, wether I can apply `lmFit` function from `limma` directly? Thanks.

Answer: Well, the easiest way is to read the data processing sections for each sample.

> The data were analyzed with Microarray Suite version 5.0 (MAS 5.0)
> using GeneData Expressionist® Pro Refiner. The trimmed mean target
> intensity of each array was arbitrarily set to 100.

The next easiest is to use [GEO's built-in analysis tools](https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc=GSE85957) (really just R scripts, but whatever) to view the value distributions. RMA normalized microarrays are typically very obvious, as their distributions are all the same.

That isn't the case here, though the distributions aren't completely nuts. These have clearly been normalized in some capacity, though you'd typically hope for bit more detail in how it was done. I don't know exactly what `limma` expects (does it expect log2 values?), so you may do a little more reading on that. You could also normalize through `limma` and just see if the results make more sense.
Question: I'm curious as to why I'm getting two different values with bedtools intersect with two sorted bedfiles:

bedtools intersect -u -a test1.bed -b test2.bed | wc -l

This gives me 157,273

However, when swapping -a and -b:

bedtools intersect -u -a test2.bed -b test1.bed | wc -l

This gives me 158,095

Answer: This is normal that order affect the result. Read carfully the help [https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html][1]

First you report all sequence from A that overlap B
And secondly you report all sequence B that overlap A

The only case where both will give you the same result it is if the two files are identical

[1]: https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
Question: I'm trying to find all the reads (by name) from a BAM file that align to various regions in a bed file. Right now I can do this with bedtools using intersectBed:

intersectBed -abam reads.bam -wo -f 1 -b regions.bed -bed

From this one can parse all the read ids that land in every interval in regions.bed, but it's not very compact. Is there a way to get bedtools to natively transform this into a more compact format, e.g.

chr1 x y .... read_id1,read_id2,read_id3

where chr1 x y is a given interval in regions.bed and the comma separated read_id1,... is the list of read ids from reads.bam that fall in that interval. In this compact format, the output BED file would have at most as many entries as there are regions in regions.bed, whereas with the -wo option it can be even larger than the number of reads in reads.bam. Thanks.

Answer: The BEDOPS `bedmap` tool with the `--echo`, `--echo-map-id`, `--delim` and `--multidelim` options can also answer this very quickly. You can pipe in reads into `bedmap` by using the `bam2bed` script:

$ bam2bed.csh < reads.bam \
| bedmap --echo --echo-map-id --delim '\t' --multidelim ',' sorted-regions.bed - \
> answer.bed

The `--echo` option returns the `region` element, while `--echo-map-id` returns a list of `read` elements which overlap the region by one or more bases.

The `answer.bed` file will look like:

[ region-A in BED format ] \t [ comma-delimited list of read IDs overlapping region-A ]
[ region-B in BED format ] \t [ comma-delimited list of read IDs overlapping region-B ]
...

As the answer is BED-formatted output, you can very easily pipe this to downstream processes that consume BED data.

Just be sure to sort regions, if not sorted, e.g.:

$ sort-bed regions.bed > sorted-regions.bed

And, likewise, for reads:

$ bam2bed.csh < reads.bam | sort-bed - > sorted-reads.bed

Or to put it together with the overall analysis:

$ bam2bed.csh < reads.bam \
| sort-bed - \
| bedmap --echo --echo-map-id --delim '\t' --multidelim ',' sorted-regions.bed - \
> answer.bed
Question: Hello all, I am new in plink analysis

I have Map file and I know that location (start /end) of ***CD2AP ***gene on chromosome 6 has base pairs 47,477,745 to 47,627,262

How can I use physical position column in MAP file to extract SNPs of this gene

Any help will be appreciated

Answer: Dear,

First you need to make sure which transcript you are interested. Keep in mind that the same SNP will have different cDNA coordinates between transcripts. Also you need to watch out for the genome version. Because g coordinates will also vary between different genome versions. Moving from your example, you can go to www.ensembl.org and write in the search CD2AP. You will see that the golden transcript is ENST00000359314 with 629 aminoacids. Now go to biomart at http://www.ensembl.org/biomart . Here select ensembl variation database, than select SNPs and Indels. Next from the left click on filters and go to gene associated filters. Here paste your gene ID which is ENSG00000198087. Next go to the attributes tab and select the coordinates you want (genomic, mRNA, coding sequence..etc) Download your results and you will have your SNPs.

I hope this helps,
Question: I'm new in this field. I want to identify circular RNAs from RNA-seq. I want to know which tools I could choose to detect circular RNAs using some simple commands. Thanks!

Answer: You could try [CIRCexplorer][1], a combined strategy to identify junction reads from back spliced exons and intron lariats using TopHat and TopHat-Fusion. The commands are very simple. For more information. please refer to [**Zhang et al., Complementary Sequence-Mediated Exon Circularization, Cell (2014), 159:134-147**][2]

[1]: https://github.com/YangLab/CIRCexplorer
[2]: http://www.sciencedirect.com/science/article/pii/S0092867414011118
Question: Hi,

I am using edgeR to carry out DE analysis for my RNA Seq samples, may I know how can I construct the design matrix and how can I make sure that during testing, the DE genes shown are of the correct contrast?

My current steps are shown below:

**Current data type:
subject: 1,2,3,4**

**treatment for each subject: T0, T1**

I would like to find DE genes between conditions (T0, T1) for each paired samples:

**design <- model.matrix(~Subject+Condition)**

**d <- d[rowSums(cpm(d)>10) >= (length(G1_ids)/2),,keep.lib.sizes = FALSE ]**

**d <- calcNormFactors(d)**

**d <- estimateGLMCommonDisp(d,design)**

**d <- estimateGLMTrendedDisp(d,design)**

**d <- estimateGLMTagwiseDisp(d,design)**

**fit <- glmQLFit(d, design, robust = TRUE)**

**result_glm <- glmQLFTest(fit)**

It seems to me that, using the above steps, I get DE genes which are differentially expressed between samples, not between conditions.

Thank you!
Maggie Pan

Answer: Please see section 3.4.2 of the edgeR vignette.

result_glm <- glmQLFTest(fit, coef=5) # Or maybe coef=4, it'll be the last column in "design"
Question: Hello,

I was wondering if there is a website in which I can retrieve the list of genes for each ontology.

For example for the biological process "Regulation of Gene Expression," I would like to have all the list that belong to this ontology etc..

Thanks in advance

Answer: I recently wrote a small tool to fetch all the children of a given term in GO: http://lindenb.github.io/jvarkit/GoUtils.html

there is an example to fetch all the associated genes:

Use GO annotation to retrieve genes associated to GO:0005216 ‘ion channel activity’

join -t $'\t' -1 1 -2 2 \
<(java -jar dist/goutils.jar -A 'GO:0005216' | cut -f 1 | sort | uniq) \
<(wget -q -O - "http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/go/gene-associations/gene_association.goa_human.gz?rev=HEAD" | gunzip -c | grep -v '^!' | cut -f3,5 | uniq | LC_ALL=C sort -t $'\t' -k2,2) |\
sort -t $'\t' -k2,2 |\
grep SCN5A -A 10 -B 10

(...)
GO:0086006	SCN2B
GO:0005244	SCN3A
GO:0005248	SCN3A
GO:0005248	SCN3A
GO:0005248	SCN3B
GO:0086006	SCN3B
GO:0005248	SCN4A
GO:0005248	SCN4A
GO:0005248	SCN4B
GO:0086006	SCN4B
GO:0005244	SCN5A
GO:0005248	SCN5A
GO:0005248	SCN5A
GO:0005248	SCN5A
GO:0005248	SCN5A
GO:0005248	SCN5A
GO:0005248	SCN5A
GO:0086006	SCN5A
GO:0086060	SCN5A
GO:0086061	SCN5A
GO:0086062	SCN5A
GO:0086063	SCN5A
GO:0005248	SCN7A
GO:0005248	SCN7A
GO:0005248	SCN7A
GO:0005248	SCN7A
GO:0005248	SCN8A
GO:0005248	SCN8A
GO:0005248	SCN9A
GO:0005248	SCN9A
GO:0005248	SCN9A
GO:0005272	SCNN1A
(...)
Question: I've ran a BLAST using the Sol genomics network database for my tomato transcriptome

http://solgenomics.net/

I have my transcriptome labeled with the sol gene IDs, but need to convert them to a more typical gene ID such as NCBI-geneID. Any thoughts on how to do this? I've tried using a series of different GeneID converters such as DAVID and biomart, but none of these programs have GeneIDs from the Sol Genomics network.

How do I convert these geneIDs? Thanks!

Answer: Okay, so I finally completed this task by doing the following:

I contacted the Boyce Thompson Institute who ran the Sol Genomics website to see if they had any conversion sheets from their Sol genomics IDs to something more common. They sent me a very good conversion sheet of Sol gene IDs to UniProt IDs.

If someone knows how to share files on biostar I'll attach the conversion sheet to my answer. It's 1.7 MB, so not too big. I can't find the same file in the Sol genomics ftp directory, but I haven't looked that hard either. so if someone knows how to share files like this or if the file is on the sol genomics ftp directory, any help would be appreciated.

So I got my sol gene IDs, got the sol to Uniprot IDs sheet, and used Galaxy to convert the sequences with the info on the conversion sheet.

Most gene enrichment analysis programs take the UniProt IDs, so this conversion did what I wanted it to do.

Also this web based program can do Gene Enrichment analysis using the Sol gene IDs: http://bioinfo.bti.cornell.edu/tool/MetGenMAP

Another good program that can use the Sol gene IDs for gene enrichment analysis: http://mapman.gabipd.org/web/guest/home

So that's how I solved my issue. Thanks for the help guys.
Question: Hi:
I have a BED file with sites of a transcription factor binding sites (mm9). I want to annotate these sites, typically located between intergenic regions, to a gene nearby. I checked a publication (PMC3080774) and it lists some tools that can map and annotate.
I tried using Chipseqanno. Initially it was straightforward till first two commands (https://bioconductor.org/packages/release/bioc/vignettes/ChIPpeakAnno/inst/doc/pipeline.html) and it went confusing.

bed <- system.file("extdata", "MACS_output.bed", package="ChIPpeakAnno")
gr1 <- toGRanges(bed, format="BED", header=FALSE)
## one can also try import from rtracklayer
gff <- system.file("extdata", "GFF_peaks.gff", package="ChIPpeakAnno")
gr2 <- toGRanges(gff, format="GFF", header=FALSE, skip=3)
## must keep the class exactly same as gr1$score, i.e., numeric.
gr2$score <- as.numeric(gr2$score)
ol <- findOverlapsOfPeaks(gr1, gr2)
## add metadata (mean of score) to the overlapping peaks

Authors did not explain why they are overlapping gr1 and gr2. Before that authors say "## one can also try import from rtracklayer". also authors did not point out names for annotations MM8 or Hg19.

My question to this forum is:
1. Given a bed file how can I map and annotate the peak regions to nearby genes.
2. Can I define 1kb, 2kb range. (apparently chipseqanno uses granges and thus is possible- however lack of mm9 and proper reference to why findOverlaps between gr1 and gr2 makes it so confusing to deal with ranges).
3. Are there any tools other than chipseqanno to achieve above two tasks.

Thanks
Adrian

Answer: You can just use *GenomicRanges* with your peaks BED file and your annotation GFF (after they have both been converted to *GenomicRanges* objects. There is a parameter in the GenomicRanges `findOverlaps()` function called `maxgap`, which would allow you to set your distances of 1kb, 2kb, etc.

With regard to why one would annotate regions that are up to 1kb, etc., away, well, that depends on your marker of study and how it is known to behave and affect gene transcription.

Kevin
Question: Hello,

I want to extract some specific strings in csv file, using either excel, R or Python.

for example as below: I want to find string from column A in column B and return in column C with 5 amino acid before and after **N**; thanks!!

```
A B C
INETTDFR MHRFLLMLLFPFSDNRPMMFFRSFIVFFFLIFFASNVSSRKQTYVIHT IVGKINETTDF
VTTSTKHIVTSLFNSLQTENINDDDFSLPEIHYIYENAMSGFSATLTDDQLDT
VKNTKGFISAYPDELLSLHTTYSHEFLGLEFGIGLWNETSLSSDVIIGLVDTG
ISPEHVSFRDTHMTPVPSRWRGSCDEGTNFSSSECNKKIIGASAFYKGYE
SIVGKINETTDFRSTRDAQGHGTHTASTAAGDIVPKANYFGQAKGLASGM
RFTSRIAAYKACWALGCASTDVIAAIDRAILDGVDVISLSLGGSSRPFYVDP
IAIAGFGAMQKNIFVSCSAGNSGPTASTVSNGAPWLMTVAASYTDRTFPAIV
RIGNRKSLVGSSLYKGKSLKNLPLAFNRTAGEESGAVFCIRDSLKRELVEGK
IVICLRGASGRTAKGEEVKRSGGAAMLLVSTEAEGEELLADPHVLPAVSLGF
SDGKTLLNYLAGAANATASVRFRGTAYGATAPMVAAFSSRGPSVAGPEIAKP
DIAAPGLNILAGWSPFSSPSLLRSDPRRVQFNIISGTSMACPHISGIAALIKSV
HGDWSPAMIKSAIMTTARITDNRNRPIGDRGAAGAESAATAFAFGAGNVDPT
RAVDPGLVYDTSTVDYLNYLCSLNYTSERILLFSGTNYTCASNAVVLSPGDLN
YPSFAVNLVNGANLKTVRYKRTVTNVGSPTCEYMVHVEEPKGVKVRVEPKVL
KFQKARERLSYTVTYDAEASRNSSSSSFGVLVWICDKYNVRSPIAVTWE
```

Answer: Using R:

getPreceding <- function(A, B, N = 4) {
x <- regexpr(A, B)
substring(B, x - N, x + attr(x, "match.length") - 1)
}

A = "INETTDFR"
B = "...IGASAFYKGYESIVGKINETTDFRSTRDAQGHGTHTAST..."
getPreceding(A, B)
# [1] "IVGKINETTDFR"
Question: I've installed BLAST and BLAT on my computer (OSX) but am confused with how many different files there are to download to set up a database on my computer and how to go about it. Can I use one with both BLAST and BLAT? Should I get the partially non-redundant database?

Answer: Once you have installed `blat`, you can download a FASTA file for your genome of interest from UCSC. If you wanted `hg38`, for instance:

$ for chr in `seq 1 22` X Y; do echo $chr; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg38/chromosomes/chr$chr.fa.gz | gunzip -c - >> hg38.fa; done

You can use this with `blat` as described in [my answer here][1].

You could run `makeblastdb` to [make BLAST databases][2]:

$ makeblastdb -in hg38.fa -parse_seqids -dbtype nucl

Or [download nucleotide databases from NCBI][3].

[1]: https://www.biostars.org/p/261839/#261983
[2]: https://www.ncbi.nlm.nih.gov/books/NBK279688/
[3]: ftp://ftp.ncbi.nlm.nih.gov/blast/documents/blastdb.html
Question: I have an assembled transcriptome that I blasted against the nr_db from NCBI. I want to use MEGAN to get some taxonomic information on the reads. I imported my blast xml files into MEGAN5 and all seemed to work fine. When I click a node to inspect which hits are associated with a given taxon I encounter a problem. The majority of the hits are listed correctly with the query id header but some do not have a header and the hole sequence or parts of a sequence is displayed. This results in incomplete sequences when I want to extract subsets of reads.

I have no clue what is happening there and how to fix it. Anyone encountered a similar problem?

Could there be a problem with the blast xml input files? Or is it a parsing problem? Any idea how to fix this?

Any help is appreciated!

Answer: I figure out what went wrong. It was a tiny formatting error in my blast input files that corrupted everything, thought blast had no problem with it. Fixed it and now everything is working fine!
Question: Can I index a BAM file that is sorted by queryname instead of coordinates using samtools?

Answer: no, indexing is only for **coordinate**