cannin/biostars_qa huggingface dataset:
Question: I have created a fasta file with 300 simulated reads at a given SNP position (reads are 75mer, and I simulate every position in the + and - strand, hence 75x4=300 reads). Then, I map them back to the ref genome using bowtie and then use pileup to see the alignment at the position in question:

```
bf <- BamFile(bamfile)
param <- ScanBamParam(which=GRanges("chr10",IRanges(start=100189138, end=100189138)))
```

`quickBamFlagSummary` shows that there are 300 reads

```
> quickBamFlagSummary(bf, param=param)
group | nb of | nb of | mean / max
of | records | unique | records per
records | in group | QNAMEs | unique QNAME
All records........................ A | 300 | 300 | 1 / 1
o template has single segment.... S | 300 | 300 | 1 / 1
o template has multiple segments. M | 0 | 0 | NA / NA
- first segment.............. F | 0 | 0 | NA / NA
- last segment............... L | 0 | 0 | NA / NA
- other segment.............. O | 0 | 0 | NA / NA

Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.
Indentation reflects this.

Details for group S:
o record is mapped.............. S1 | 300 | 300 | 1 / 1
- primary alignment......... S2 | 300 | 300 | 1 / 1
- secondary alignment....... S3 | 0 | 0 | NA / NA
o record is unmapped............ S4 | 0 | 0 | NA / NA
```

Pileup only shows 262 reads (63+62+75+62)... Any idea why? I don't use any threshold on `min_base_quality` and `min_mapq`...

```
> pileup(bf, scanBamParam=param, min_base_quality=0,min_mapq=0)
seqnames pos strand nucleotide count which_label
1 chr10 100189138 + A 63 chr10:100189138-100189138
2 chr10 100189138 - A 62 chr10:100189138-100189138
3 chr10 100189138 + G 75 chr10:100189138-100189138
4 chr10 100189138 - G 62 chr10:100189138-100189138
```

I looked the region up with IGV, and also, IGV reports 300 reads.

I actually simulate several different positions, and the 262 reads (instead of 300) appear every time.. Not sure why. I am sure pileup is filtering the reads somehow..

Also, the same exact script works fine if I simulate 36mer reads (instead of 75...)! Any ideas why?

Answer: Thank you! it is solved by changing `max_depth` parameter within `PileupParam`()

```
> pileup(bf, scanBamParam=param, pileupParam=PileupParam(max_depth=8000, min_mapq=0, min_base_quality=0))
seqnames pos strand nucleotide count which_label
1 chr10 100189138 + A 75 chr10:100189138-100189138
2 chr10 100189138 - A 75 chr10:100189138-100189138
3 chr10 100189138 + G 75 chr10:100189138-100189138
4 chr10 100189138 - G 75 chr10:100189138-100189138
```

I was adding my arguments without PileupParam constructor, this is wrong:

```
> pileup(bf, scanBamParam=param, min_base_quality=0,min_mapq=0, max_depth=8000)
seqnames pos strand nucleotide count which_label
1 chr10 100189138 + A 63 chr10:100189138-100189138
2 chr10 100189138 - A 62 chr10:100189138-100189138
3 chr10 100189138 + G 75 chr10:100189138-100189138
4 chr10 100189138 - G 62 chr10:100189138-100189138
```
Question: Dear all,

This is a newbie question :)

I'm building a linear model to identify significant predictors of mutation count/types in tumours from TCGA. I want to include expression levels of a couple of genes, but I am quite new to RNA-Seq analyses and best practices. TCGA provides RNA-Seq data at the gene level in three formats: HTSeq-counts, FPKM and FPKM-UQ. I have been reading (tutorials and the questions here) and asking around and I have reached the conclusion that I can use FPKM-UQ values to compare across samples without any further pre-processing - Is this true? Or would you recommend doing pre-processing to these values before comparing?

Thanks so much,
Daniela

Answer: Hi Daniela,

I would highly recommend the HTSeq counts, actually, because these will be raw counts. The normalisation method that produces FPKM expression levels has come under criticism in recent years and is now not even recommended by some sources. The main issue with this [FPKM] method is that cross-sample normalisation is non-existent, as such, it's akin to comparing multiple batches without even doing any correction for batch.

Use HTSeq counts and load these into DESeq2 or EdgeR for downstream analyses.

I have recently analysed an entire TCGA RNAseq dataset (>500 samples) and I used HTSeq counts. They work very well.

Kevin

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

**Update May 2, 2018:**

The TCGA states that "*To facilitate cross-sample comparison and differential expression analysis, the GDC also provides Upper Quartile normalized FPKM (UQ-FPKM) values and raw mapping count.*" - https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification

My original advice still stands, i.e., better to obtain the raw HT-seq counts (where available), and re-process those using an updated normalisation method, like TMM (EdgeR) or geometric (DESeq2). Some TCGA datasets are only available in RSEM counts, which are also possible to use and input to DESeq2 using *tximport*
Question: I am looking through the Biopython Tutorial and Cookbook (http://biopython.org/DIST/docs/tutorial/Tutorial.html) and in order to look at DSSP files with biopython, I need to obtain a copy and license of DSSP from this link (https://swift.cmbi.umcn.nl/gv/dssp/). Could somebody please explain how I take the files from github and install them to my computer? I downloaded the files on my computer and put them in my current working directory in python and I'm still not able to use DSSP files in biopython

Answer: I'd suggest using `conda` to install it:

`conda install -c salilab dssp`

So long as you have the `dssp`* binary available in your `$PATH`, BioPython should be happy.

*It may be called `mkdssp`
Question: I need to design a primer set for a gene from cDNA. How can I go about making a primer set if it doesn't matter what portion of the gene I use? In addition, how can I predict the size of the application product resulting from using the primers I design?

Answer: I'll assume:

- You know about making primers for gDNA targets
- You are designing a primer for RT-qPCR (since you didn't explicitly state that)
- My suggestions below mainly are based on experience for human genes but are most likely transferable to other eukaryotes (but others may want to add specific things about other species).

I would suggest doing the following:

1. Get the mRNA sequence of your gene of interest, from UCSC or Ensembl, doesn't matter
2. Use this sequence in [Primer3][1]
3. The best is to have either a primer on an exon-exon junction, e.g. primer partially in exon 2 and exon3. This will limit the chance of amplification of contaminant gDNA. In Primer3 you can mark this junction site using '-' in the sequence or use the box "Overlap Junction List". The explanation on the website is quite clear.
4. Second best would be an amplicon that spans a (large) intron, so one primer in exon 4 and the other in exon 5. This again will reduce the chance of having contaminant gDNA amplification. In this case, you would make the junction a target. In Primer3 you would mark the sequence of your target with '[' and ']' or use the "Targets" box.
5. Test the primer pair using [in silico PCR][2] and [blat][3], both on the genome assembly and the Gencode genes to investigate potential off-target hits. If any hits are found (and their amplicons are quite small) repeat the process to get other primers, for example using Masking of the primer site which didn't turn out to be a good one. In Primer3, masking can be done by changing the sequence to N, or excluding the sequencing using '<' or '>', or use the "Excluded Regions" box.

[1]: http://primer3.ut.ee
[2]: https://genome.ucsc.edu/cgi-bin/hgPcr
[3]: https://genome.ucsc.edu/cgi-bin/hgBlat
Question: Hi,
I am currently getting the following error:

library(biomaRt)

mart<-useMart(biomart="ENSEMBL_MART_ENSEMBL",
dataset="hsapiens_gene_ensembl",
host="grch37.ensembl.org")

Extra content at the end of the document
Error: 1: Extra content at the end of the document

I can't find any post specifying any scheduled maintenance. Is the server down?

Thanks!

Answer: Hello,

There was a temporary disruption with the servers at the time you experienced this error message yesterday. Everything is running as it should be now. If you encounter any further errors (or are still encountering the same problem), please e-mail helpdesk@ensembl.org, and we will be able to assist you from there.

Best wishes

Ben
Ensembl Helpdesk
Question: Hi,

Given,I have a gene expression dataset. I want to find out, all the genes that are highly correlated with the given gene/ have the same expression pattern.

Example, if there would be an R package, it should take one input as expression data matrix and another input as gene of interest and should provide an output of a list of genes, that could have same expression pattern or are highly correlated.

I know that I could write a function in R regarding this and filter the genes on the basis of their correlation values and take into account the cor value more than 0.5. But, it takes so much of time, when you have a data with more than 20,000 rows and 1,000 columns.

Answer: Using the function below is pretty fast using 1500 columns and 30k rows

```
mat <- replicate(1500, rnorm(30000))
gene <- mat[sample(1:30000,1),]

system.time({
resp <- c()
est <- c()
for (i in 1:30000){
x <- cor.test(gene, mat[i,])
if (abs(x$estimate[[1]]) > 0.5){
resp <- c(resp, x$p.value)
est <- c(est, x$estimate[[1]])
}
}
})

user system elapsed
7.69 0.02 7.71
```
Question: I've got a few thousand small bam files produced against the exact same reference, and I want to merge them into one single big bam file. What is the best way to do that?

Should I do this iteratively or can I pass a long list of bam files to samtools/picard/etc in one go?

Edited, since this is now partially solved. In my terminal, the methods below works for up to 4092 files. More than that raises an error:

samtools merge all.bam *.bam
samtools merge all.bam `find /basedir/ -name "*myfiles*.bam"`
samtools merge all.bam /basedir/*/???/*myfiles*.bam

Answer: Bamtools handles this very cleanly. Note that this properly constructs the header, which might matter if you have different read groups in each file.

bamtools merge -list files.bamlist -out merged.bam

You can also supply a region via the -region flag, using the format 12:2232..3328.
Question: I would like to download the list all SNPs for a gene of my interest with all their details like Chromosome, Physical position etc. How can I do this?

Answer: Dear Mahantesh,

For people who are not used to command line, biomart has a very nice interface for downloading SNPs etc. Go to http://www.ensembl.org/biomart/martview/ . You can filter based on gene IDs, you will find these on ensembl/UCSC.
Question: H!
Which bioinformatics tools supporting CRISPR-Cas Experiments do you recommend? Or which tools are you using?

Answer: Hi
You can try our package in R, [Casilico][1]. It supports Cas12, Cas13, and Cas14 CRISPR-Cas systems. Using a list of important features such as mismatch tolerance rules, self-complementarity, GC content, frequency of cleaving base around the target site, target accessibility and PFS (protospacer flanking site) or PAM (protospacer adjacent motif) requirement, CaSilico searches all potential crRNAs in a user-input sequence.

CaSilico outperforms previous crRNA design tools in the following respects: 1) supporting any reference genome/transcriptome for which a FASTA file is available; 2) designing crRNAs that simultaneously target multiple sequences through conserved region detection among a set of sequences; 3) considering new CRISPR-Cas subtypes; 4) reporting a list of different features for each candidate crRNA, which can help the user to select the best one. Our paper was puplished in [Frontiers in Bioengineering and Biotechnology][2]. ![enter image description here][3]

Our tool is sufficiently flexible to tune some key parameters governing the design of crRNA and identification of off-targets, which can be led to increases the chances of successful CRISPR-Cas experiments. Given these capabilities, Casilico addresses end-user concerns arising from the use of sophisticated bioinformatics algorithms and has a wide range of potential research applications in different areas especially design of crRNA for pathogen diagnosis. Casilico was successfully applied to design crRNAs for different genes in SARS-CoV-2 genome, as some of the crRNAs have been experimentally tested in the previous studies.

[1]: https://github.com/mrb20045/CaSilico
[2]: https://www.frontiersin.org/articles/10.3389/fbioe.2022.957131/full
[3]: https://user-images.githubusercontent.com/9910942/158050606-9c592f1c-c0a4-4a4e-8cf8-ded8e0c7e7b6.png
Question: Hello,

I took a 10x matrix from a collaborator and created a Seurat object. The issue I am having which I only realised when attempting to visualise my favourite genes, is that the original matrix has gene names in the format of "gene name - ensembl"

example :

Dcn-ENSMUSG00000019929.16, Inmt-ENSMUSG00000003477.5, Mfap4-ENSMUSG00000042436.12

So I can not search for my genes of interest unless I know what the corresponding ENSMUS following it is.
Is there I way I can remove at this stage everything after (and including) the "-" so that I am left only with the gene name?

So convert "Dcn-ENSMUSG00000019929.16" to simply "Dcn".

Thank you :)

Answer: Huh, you are quite right that my comment on jomo's answer doesn't work. Too used to `SingleCellExperiment` objects, where that does work. Anyway, there are two workarounds. The annoying one is to extract the counts from assay data and remake the Seurat object:

# Get counts
count.data <- GetAssayData(object = pbmc_small[["RNA"]], slot = "counts")
count.data <- as.matrix(count.data)

# Rename rows
rownames(count.data) <- unlist(map(strsplit(rownames(count.data),"-ENSMUS",fixed=T),1))

# Generate new Seurat object.
new.obj <- CreateSeuratObject(
count.data,
project = "SeuratProject",
assay = "RNA",
min.cells = 0,
min.features = 0,
names.field = 1,
names.delim = "_",
meta.data = NULL
)

Obviously not ideal. You'll have to redo the normalization. Other option is to convert to a `SingleCellExperiment` object, rename the rownames as above, then convert back to a Seurat object. This at least retains the normalization and any dimensionality reductions, metadata, clustering, etc.
Question: I have a table:

W X Y Z
A 2 3 4
A 2 3 6
B 1 2 3
C 3 2 1
B 1 3 4
B 1 2 2

# reproducible data
table <- read.table(textConnection("W\tX\tY\tZ
A\t2\t3\t4
A\t2\t3\t6
B\t1\t2\t3
C\t3\t2\t1
B\t1\t3\t4
B\t1\t2\t2"), header = TRUE)

I want to combine all rows with same value in Column W, so that it looks like

W X Y Z
A 2 3 2,6
B 1 2,3 3,4,2
C 3 2 1

**Note:** All different values of column **W** entries appear as comma separated.

Answer: You can also use the `aggregate` function in combination with the `unique` function

aggregate(table[,2:4], list(table[,1]), function(x) paste0(unique(x)))
Question: can someone point to a detailed explanation of how UCSC genome browser interprets coordinates for, and displays, minus strand sequences? It's very confusing. For example this genomic GFF coordinate from mm9 genome:

chr1:3206080-3206102:-

Pulls out the sequence:

>seq1
AATACAAGGAGCCGCAATGTCCA

which is correct. However, in the genome browser it appears as the reverse complement going from right to left

>seq2
TGGACATTGCGGCTCCTTGTATT

Questions are:

why does it UCSC reverse complement it and display it right to left in this way? What is the logic behind this? Is there a way to reverse it, such that the sequence displayed is the sequence pulled from the genome (i.e. seq1)

I thought UCSC is always BED coordinate based, which is a 0-based coordinate system not 1-based. In that case I would have expected the sequence determined by chr1:3206080-3206102 (meant to be a GFF coordinate) to be one base shorter than the sequence I get from UCSC, i.e. chr1:3206080-3206102:- in BED should have a minus 1 start, yielding chr1:3206079-3206102:-

any intuitive explanation of how to think about UCSC orientation display choices and how that connects to sequences and their reverse complement on minus strand (and where/when BED versus GFF conventions are used) would be very helpful.

Answer: Huh. Yeah, that's not well documented, you are right. I looked in a few places where I expected to see it, but can't find a place to give you a handy link on the viewer reverse strategies. I'll tell you what I know, but that is not the official word and they may have better details about the underlying decisions. And it may just be that I haven't had enough coffee to find it in the documentation this morning.

On the graphical display: Yes, it displays the "upper" strand 5'-3' by default. There are 2 ways to flip it. One is an old way and slightly hard to find. The other is the "reverse" button on the mid-page buttons. That's newer and better, I think, because I really like the way it load the track labels to the other side to remind you that you have flipped the graphical view. But I'll show you both on the image below.

Old way: when you turn on the "base position" track to "full", and you are zoomed in enough to see the 3 frame translation, you will also have access to a little teeny arrow. Clicking that arrow will flip the display.

New way: click that mid-page "reverse" button. The flipped view isn't shown below, but you'll see it immediately if you just click that.

On the 0-1 based coordinates, welcome to the club, you are now officially one of us. This bites everyone in the a** at some point in their bioinformatics career. Here's their documentation on that:
http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1

So part of the answer is just: that's the way it is. But there are tricks. Does that help?
Question: Hi guys,

I've got a question : what's the main difference between the mapping and the quasi-mapping? (I used Salmon)

Actually, I read many papers about the subject but I'm still struggling with that.

Both of these kinds of mapping use a suffix array and a hash table (k-mers) to increase the mapping speed.

I also know that quasi-mapping (like Salmon) kills two birds with one stone by mapping and couting reads for each transcripts

So, where is the difference if mapping and quasi mapping use a suffix array and a hash table

Best,
Vincent

Answer: Hi @vincentpailler,

The main difference is in how much computational effort is expended on identifying the optimal alignment of each read. Basically, for certain tasks (e.g., variant detection), it's crucial to know for each read, both where it arises from in the underlying reference, and how it aligns to that reference sequence (since, e.g., repeated mismatches will be used as evidence of the existence of a variant). However, for many problems, including RNA-seq quantification, that precise information isn't necessary. If we know which transcripts a read arises from, and the likely position and orientation of the read on that transcript, then we can use that information for accurate abundance estimation. So, the main differences are in how much effort is expended to find the optimal alignment between the reads and the reference. Additionally, when aligning to the transcriptome (rather than the genome), there are a lot of repetitive mappings, since alternative splicing gives rise to many identical / near-identical alignments for a read. That is, the read may align to different transcripts and positions, but these are produces as a result of alternative splicing, and so the underlying reference sequence is identical. This leads to a lot of redundant work solving identical alignment problems repeatedly. Salmon's mapping algorithm is optimized to deal with such repetition. Recent version even include a `--validateMappings` flag which (optionally) _does_ perform an alignment of the read at the candidate location --- however, to keep this process efficient, it is necessary to track when reads map to identical underlying reference sequence so as to avoid redundant work.
Question: Hi,

in IGB, it is possible to access directly a lot of data, included ENCODE datasets (from the Data Access tab). Is it possible to search for a dataset by its name (instead of browsing the while tree)?

thanks,

Arnaud

Answer: Hello Arnaud,

My name is Mason Meyer, support specialist for the IGB project. Thank you for taking the time to post this excellent question.

Currently, it is not possible to search the tree, but this is actually a feature planned for IGB 9.0.0, our next major release. In the meantime, we would love to hear any additional suggestions you may have for this feature, or even other features you would like to see make it into IGB 9.0.0.

Thanks again Arnaud!
Question: I have a large WGS dataset composed of ~350 individuals from 6 different populations and I want to figure out the fixed differences between each population (and the rest) for the sake of calculating the direction of selection: DoS = Dn/(Dn + Ds) - Pn/(Pn + Ps).

I can easily calculate the within-population polymorphic differences (Pn/Ps) by splitting the VCF based on population and running each through SnpEff separately, but it is the Dn/Ds ration I'm having trouble figuring out (fixed differences between a given population and the rest). Any ideas?

Answer: `bcftools` will do this with `--private`. See [here][1]. You'll need to specify the samples belonging to the population of interest.

[1]: https://samtools.github.io/bcftools/bcftools.html
Question: For my RNA-seq project, I have a list of genes in one column and I am trying to associate this list with a second column that has all the biological pathways that each respective gene is involved in. I looked at databases such as PathwaysCommons, but they don't seem to provide a non-network based representation that I can FTP directly into my local computer.

Is there a way to get a list of genes and their respective associated biological pathways? One gene can have many very many pathways that it is known to be involved in. Is there any online utility that I can use to find this kind of information?

Answer: I suggest that you try and do this by using either:

- [BioMart][1]
- [DAVID][2]

They offer a GUI that allow you to easily append pathway information to a list of genes.

[1]: http://www.ensembl.org/biomart/
[2]: http://david.abcc.ncifcrf.gov/
Question: Hi all,

I'm about to align a set of short reads on the genome after an exome capture.

I wonder if it would be worth to replace all the non-genic genomic sequences with 'N'.

Wouldn't it speed up the alignment with BWA ? What would be the other possible consequences ? Has it ever been done before ?

Many thanks

Pierre

Answer: In a word, no. It's not necessary for any of the BWT aligners. In fact it would cause you to miss any unannotated exons, or possibly entire novel genes. It would also hide other biological events, such as genomic sequence contamination in the preparation process.

If you want to speed up the alignment, split the reads into several chunks and map them in parallel. Afterwards use a merging program such as samtools or Picard to create a single BAM file. We usually aim for ~ 1Gbp per process for a vertebrate genome.
Question: Hi All,

can we create abstract image or pattern specific to person's DNA sequence ?

design something similar to : http://bit.ly/1SjCewK

Thanks,
Jay

Answer: Yes, I would recommend using a WebGL framework like [three.js][1] like in this [example][2] where they visualize music or this [example][3] for a more "realistic" look.

[1]: http://threejs.org/
[2]: http://macobo.github.io/WebGL-Audio-Visualization/#/particles
[3]: http://oos.moxiecode.com/js_webgl/dna/