cannin/biostars_qa huggingface dataset:
Question: Hi,

For running [ActivedriverWGS](https://github.com/reimandlab/ActiveDriverWGS) software I will need coding or non coding parts of genome in BED12 format. I have found coding part of genome (in txt format though). But I don't know how to find non coding part of genome (BED12 format) also Transcription factor binding in BED4 format. I have contacted the developer but no response. Any suggestion please?

Answer: Sorry, finally I got what I want

I actually need enhancers, promoters, or other regulatory elements of human genome to find cancer driver genes placed in these regions

People says by

awk '$3=="transcript" &&
$20!="\"protein_coding\";" &&
$20!="\"translated_processed_pseudogene\";"' gencode.gtf

Will return non-coding parts of regions

like

awk '$3=="transcript" && $20!="\"protein_coding\";"{print $20}' gencode.gtf | sort | uniq -c | sort -nk1
1 "translated_processed_pseudogene";
2 "Mt_rRNA";
3 "IG_J_pseudogene";
3 "TR_D_gene";
4 "TR_J_pseudogene";
5 "TR_C_gene";
10 "IG_C_pseudogene";
18 "IG_C_gene";
18 "IG_J_gene";
22 "Mt_tRNA";
25 "3prime_overlapping_ncrna";
27 "TR_V_pseudogene";
37 "IG_D_gene";
58 "non_stop_decay";
59 "polymorphic_pseudogene";
74 "TR_J_gene";
97 "TR_V_gene";
144 "IG_V_gene";
182 "unitary_pseudogene";
196 "IG_V_pseudogene";
330 "sense_overlapping";
387 "pseudogene";
442 "transcribed_processed_pseudogene";
531 "rRNA";
802 "sense_intronic";
860 "transcribed_unprocessed_pseudogene";
1529 "snoRNA";
1923 "snRNA";
2050 "misc_RNA";
2549 "unprocessed_pseudogene";
3116 "miRNA";
9710 "antisense";
10623 "processed_pseudogene";
11780 "lincRNA";
13052 "nonsense_mediated_decay";
25955 "retained_intron";
28082 "processed_transcript";

But I am not sure from these regions which parts are related to `enhancers`, `promoters`, or `regulatory elements`
Question: Dear All,

It might be a naive question, but I got stuck a bit. How could I read csv file into R and make the gene column had no label like the picture in attached. everytime I read my data in, the label is still there.

thank in advance.

regards,

synat

![enter image description here][1]

[1]: /media/images/f7f0f656-9dc3-431b-992f-c9ed9b85

Answer: I'm going to interpret your question a bit. I'm assuming here that you just want to convert your `data.frame` to a matrix.

```
# create a matrix from the data.frame
mat2 = as.matrix(mat[,-1])

# add rownames to the matrix
rownames(mat2) = mat[,1]

# show the matrix with genes as rownames
mat2
```

Note that if the "genes" column has duplicates, this won't work as expected since rownames must be unique.
Question: Hello there,

I am having a lot of MiSeq data I am trying to analyze and I figured out using FastQC that I have a lot of fails in the report it generates and I wanted to pick your brain to get sense of what should be done in that case.

As you all know FastQC generates this kind of information:

```
[PASS] Basic Statistics
[FAIL] Per base sequence quality
[PASS] Per sequence quality scores
[FAIL] Per base sequence content
[FAIL] Per base GC content
[WARNING] Per sequence GC content
[PASS] Per base N content
[WARNING] Sequence Length Distribution
[FAIL] Sequence Duplication Levels
[WARNING] Overrepresented sequences
[FAIL] Kmer Content
```

The issue here is that I am analyzing targeted sequencing data, so I am expecting to have a lot of duplications, what I don't clearly see is whether or not to take the result of FastQC as correct based on the standard they are publishing on their website (how a good report should look like and how a bad one should look like), so I am expecting the GC content to go crazy with the amount of duplication because of the type of experiment (deep sequencing), now based on the information provided in the example above, do you think fastq post processing like clipping and trimming would correct the reads or is it failing in the level of the MiSeq machine already (experimental contamination?)

Rad

Answer: I consider FastQC as a screening program for rapidly getting QC info; I don't see it as a definitive QC analysis tool, and it's not designed to be that. It pulls a subset (10% IIRC) of the reads in the file and works with those. It's designed to be very picky and give false positives. With a yellow or red flag, FastQC is basically saying "Yo! Look at this! Might be a problem in a general context." If you see that the result is in fact what you expect given the experimental context, then it's probably safe to proceed.
Question: Hi,

I am currently analysing a RNAseq dataset of human lymphocytes containing a vector (which has mCHERRY as reporter gene and a different gene A).

In order to include the annotation for the new genes I have followed the instructions of Alexander Dobin (https://groups.google.com/forum/#!topic/rna-star/FGQRotrCB1Q) and it seems like STAR (v2.7.3a) is detecting them and including them in the BAM files.

However, when I try summarising the BAM files using featureCounts (I am using Rsubread 1.34.4 in my university's computer cluster and subread 2.0.0 at my local machine) these genes are not listed (not even with 0 counts). I have also tried adding an additional line for each gene changing "exon" for "gene", but it does not solve the issue.

Do you have any suggestions of why this might be?

This is the code I am using in the Rsubread:

feature_counts <- featureCounts( nthreads = 6,
isGTFAnnotationFile = TRUE,
GTF.featureType = "exon",
annot.ext = "genome/gencode.v31.primary_assembly.annotation.gtf",
GTF.attrType = "gene_id",
primaryOnly = TRUE,
isPairedEnd = TRUE,
files = paste( "alignment/",
list.files( path = "alignment/",
pattern = ".bam$",
recursive = TRUE,
include.dirs = TRUE )[-1],
sep = "" ) )
And in subread:

featureCounts -T 2 -t exon -a ./gencode.v31.primary_assembly.annotation.gtf -t exon -g gene_id -M --primary -p -o ./count_matrix/counts.txt \./alignment_test/AG*/*.bam

Many thanks!!

Eduardo

Answer: After checking a bit more carefully at the GTF file I found that even though in the instructions from Alenxander Dobin it suggests adding a space between the start and end of the gene (e.g. 1 and 708) and this seems to work for STAR, this should be a tabulation and not a space. After changing this it seems like featureCounts has no problems with including the new genes in the summary.
Question: Many experts claim that the negative binomial distribution is better than the Poisson distribution for modeling discrete RNA-Seq data.

Is this because

The expression profiles of organisms resemble the negative binomial distribution - i.e.if I binned real gene transcription counts from an RNA-Seq experiment it would look like the plot below.

or

The sampling error of gene expression is such that the true population mean of a gene looks like the negative binomial - i.e the true mean expression level of that gene is probably more (because of the skew) than the mean expression level of a sample of reads of that gene drawn from replicates.

or are these two concepts the same thing?

Answer: It's not about the skew of the distribution, it's about the variance. Take a look at the top two panels of this image (self link). The big difference is that the distribution based on a NB distribution can be shorter and fatter, or overdispersed.

The key to understanding the negative binomial distribution is that it's the same as taking lots of poisson distributions with slightly different means and adding them together. Essentially, this is saying that there's some sort of bias in the data causing this variability, but we aren't sure exactly what's causing it. The obvious candidates (GC content, mapability, sequence composition, etc) don't explain it all. Since it's difficult to correct for something when we don't know the source, we can do the next best thing, which is to model the distribution we see accurately. That's where NB distributions can overcome some of the limitations of poisson (namely that the variance/mean ratio is always equal to one).
Question: Dear all,

I am running this command:

Rsubread::buildindex(basename='output/genome/ref/index/subread/celegans',reference="output/genome/ref/seq/celegans.fa")

to build an index for the C. elegans genome. R then outputs lots of interesting information to the console.

Do you know how I can redirect that output to a file?

Thanks.

C.

Answer: sink("yourOutput.txt")
your stuff
sink()
Question: Hi all,
I'm playing with an alignment algorithm, that should look like BLAT (I think).

From the Reference sequence, I've sorted all the words of size=K

AAAAAAACT index in ref=1901
AAAAAAACT 2504
AAAAAACTT 2505
(...)
TTTTTTACC 189
TTTTTTTTA 87

when we want to map a sequence on the chromosome, the query is also indexed

AAAAAAACG 20
AAAAAACGA 21
(...)
TTAATGGAG 4

so, the words of the query can be quickly mapped to the reference using a binary-search.
At the end, I got an array of hits:

struct Hit
{
int position_in_ref;
int position_in_query;
int word=K;
}

that should look like this (the figure is from the BLAT paper)

Now, how should I use this information to get the best alignment(s). I suspect a simple solution should be something like a dynamic algorithm but I'm lost at this point. Any suggestion ?

Answer: So you've done your search stage and obtained candidate alignment seeds. Now you need to choose groups to join and fill the gaps. How to select groups of seeds? BLAT merges overlapping seeds to start with. Then you could look for groups that form chains within a certain diagonal band. You can find the band by subtracting the query base coordinate from the reference base coordinate. Then sort everything by band and the chains will be grouped together. I think BLAT recurses at some point, to further join these seeds by repeating the Kmer finding operation locally, with smaller and smaller Kmers.

Or, if you don't want to recreate BLAT exactly, you can use dynamic programming (e.g. Smith Waterman) to stitch them together. If your seeds are good enough, you could assume they are true and use DP to join them up, or you could go for it and DP re-align the whole area to get the best possible result, using banded alignment. A banded alignment means that you can ignore the DP cells outside your diagonal band and speed up the calculation. When you say "best" alignment, that implies to me a rigorous search i.e. SW, but maybe I've misunderstood that point.
Question: Hi,

I am trying to analyze a GSE that was performed on affymetrix HT HG-U133+ PM Array Plate [GPL13158][1]

But I cannot find the corresponding annotation package on bioconductor.

Affymetrix website is of no use, return a 404 error when selecting the array.

If anyone has experience with this array and know how to map gene name to probe name, I would be glad to have some help

Thanks

[1]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL13158

Answer: Oh yeh, some of those links to Affymetrix pages will be broken due to [I believe] the merger between ThermoFisher and Affymetrix. You can now find technical information on either or both websites, like, for example: HT HG-U133+ PM Array Plate - Support Materials.

The *HT_HG-U133_Plus_PM* is merely a modified version of the original U133 - it has no mismatch (MM) probes and has a few other modifications.

There is a Bioconductor package HERE. You may have to install it from the source tar.gz.

If not the Bioconductor package, then just download the annotation table at the bottom of the GEO page (GPL13158), and load that manually into the fdata slot of the *ExpressionSet* object.

Kevin
Question: I am using `edgeR package` in `R` and here the code I used

library(edgeR)
target <-read.delim("metadata.txt", header=T)
x <- read.csv('data4.txt', row.names = 1)
group <- factor(target$group)
design <- model.matrix(~0+group)

when I print the design variable I would get the following results:

> design
group1 group2
1 1 0
2 1 0
3 0 1
4 1 0
5 0 1
6 1 0
7 1 0
8 1 0
9 0 1
10 0 1
11 1 0
12 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

I would like to make a new variable called comparison like this:

comparison = "group1-group2"

I actually do not know how to do that because group1 and group2 are not headers or separated sections. do you guys know how to that?

Answer: Get the dimnames:

#example data
group <- factor(1:2)
design <- model.matrix(~0+group)
design
# group1 group2
# 1 1 0
# 2 0 1
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$group
# [1] "contr.treatment"

paste(dimnames(design)[[2]], collapse = "-")
# [1] "group1-group2"
Question: Hi there!

I got a list of proteins ([UniProt ][1]ID) from a proteomics experiment and i am trying to build a model containing all the routes where my proteins are implicated. To achieve this, my first idea is to extract the [Reactome][2] route identifier from all the Uniprot IDs and then obtain the SBML for each one of those routes.

Joining the models (SBMLs) is a problem for future me, but i actually dont know how to obtain the SBML from [Reactome][2] with a list of identifiers like this (R-HSA-112409). I obtained the list of identifiers with a small python script.

I'm trying to build the model, not to get the routes more represented within my list of proteins. [Reactome][2] provides this type of analysis, but the SBML resulting from the top routes does not contain the same information that those accesed from the UniProt identifiers.

This is my first contact with SBMLs, so i would like to ask for advice for constructing the model (other databases providing SBMLs given a list of proteins for example), and if it were possible, for a way to obtain the SBMLs from [Reactome][2] using the route identifiers!

Thank you for your answers!

[1]: http://www.uniprot.org/
[2]: http://www.reactome.org/

Answer: I think you can get what you want with the [REST API][1]. For example, for R-HSA-112409:

curl http://reactome.org/ReactomeRESTfulAPI/RESTfulWS/sbmlExporter/112409 > R-HSA-112409.sbml

[1]: http://reactomews.oicr.on.ca:8080/ReactomeRESTfulAPI/ReactomeRESTFulAPI.html
Question: **Can I transfer this first format to the second one just by basic shell procession or awk or sed on linux? This is a toy example:**

This kind of text file is what I have, three cols, col2 and col3 like range, left close and right open,
chr1 0 2 0
chr1 2 6 1.5
chr2 0 3 0
chr2 3 10 2.1

Transfer to describe each position as:
chr1 0 0
chr1 1 0
chr1 2 1.5
chr1 3 1.5
chr1 4 1.5
chr1 5 1.5
chr2 0 0
chr2 1 0
chr2 2 0
chr2 3 2.1
...
chr2 9 2.1

Someone has idea how to solve? Thanks!!!

Answer: If you want to use set operations:

$ bedops --chop 1 foo.bed | bedmap --faster --echo --echo-map-id --delim "\t" - foo.bed | cut -f1,2,4 > answer.txt
Question: Hi Guys,

I have this table (mydf) with the column ALT, and I want to split the string in each cell of this column so that there are only one character per cell and `cbind` them to new column as ALT1, ALT2, ALT3, ..., ALTn based on the number of characters in the cells and get the result table. How can I do this in R? Thank you for your help.

`mydf`

```
REF ALT GENE
A GC MAPK
T GCA..n MAP2K
```

Result:

```
REF ALT ALT1 ALT2 ALTn GENE
A GC G C MAPK
T GCA..n G C nth character MAP2K
```

Answer: I think I figured it out myself:

```r
library(plyr)
sst <- strsplit(mydf[,"ALT"], "")
additional.cols<-rbind.fill.matrix(lapply(sst, rbind))
colnames(additional.cols)<-c(paste("ALT",sequence(ncol(additional.cols)),sep=""))
```
Question: Hi all,

I'm new with using MACS and HOMER.
I have created a bed file using MACS, and got 5 columns, for example:

chr1	770811	771252	MACS_peak_1	93.92
chr1	851512	851790	MACS_peak_2	57.31
chr1	920907	921486	MACS_peak_3	207.37
chr1	946784	946982	MACS_peak_4	73.24

Does anyone knows how can I create a column that contains the Unique Peak ID for HOMER analysis?
Also, What is the difference between using `findMotifs.pl` and `findMotifsGenome.pl` ? besides the fact that the first uses Unique Peak ID and the second uses genomic positions?

Appreciate your help!

Answer: The difference between `findMotifs.pl` and `findMotifsGenome.pl` is that the first works with gene names (or IDs) as input, and the second with genome coordinates. You have coordinates so you can use `findMotifsGenome.pl`.

For this you need a bed file, I think only the first 4 columns are needed (I am not sure though if the 5th is ignored automatically, or that you have to remove it yourself first). The unique peak ID is already in the 4th column.

Something like this will find the motifs of your peaks.

findMotifsGenome.pl File_peaks.txt ~/Reference/hg38/ensembl/Homo_sapiens.GRCh38.dna.primary_assembly.fa File_peaks -p 20
Question: Dear Users,

I couldn't find any script that could split my vcf file that has scaffolds (10322) but not chromes number, i would like to split each scaffold into individual vcf file or a text file.

Thank you

Answer: First get a list of uniq scaffolds:

bgzip merge-test.vcf
tabix merge-test.vcf.gz
gzcat merge-test.vcf.gz | grep -v "^#" | cut -f1 | sort | uniq > scaff_names​

Then Split by individual scaffolds:

parallel -a scaff_names "bcftools view merge-test.vcf.gz {} > {.}.vcf"​

or by for loop:

for scaff in `cat scaff_names`; do bcftools view merge-test.vcf.gz $scaff > $scaff.vcf; done​
Question: Hello everyone,

I have a simple question, the file `series_matrix.txt` in GEO are always the normalized data of the study?

Best,

Leite

Answer: Hey Leite,

The answer is that, yes, the series matrix files should contain normalised, log2 values. However, the GEO provide situations in which these files may not contain normalised data:

> GEO2R operates on **Series
> Matrix files which contain data extracted directly from the VALUE
> column of Sample tables. Submitters are asked to supply normalized
> data in the VALUE column,** rendering the Samples cross-comparable.
> **The majority of GEO data do conform to this rule.** **GEO applies no further processing other than to perform a log2 transformation on**
> **values determined not to be in log space** (see Options section).
> However, some studies, such as dual channel loop design data, may
> generate values that do not have a common reference and are not
> directly comparable. Some studies may contain Sample value data that
> are not normalized, or have a design such that the Samples were never
> intended to be directly compared. Yet other studies do not have
> sufficient replicate Samples to perform a robust statistical analysis.
> **Users should examine the original Series to understand the**
> **experimental design, and check the 'Data processing' field or VALUE**
> **description in the original Sample records for information on what the**
> **values represent**. **The box plot feature on the Value distribution tab**
> **is provided to help users assess whether the distributions of values**
> **across Samples are median-centered, which is generally indicative that**
> **the data are normalized and cross-comparable.**

*[source: https://www.ncbi.nlm.nih.gov/geo/info/geo2r.html]*

When you obtain data, you should always check the distribution with box- and scatter plots, and histograms, in order to gauge whether thy are normalsed or not.

Kevin
Question: Hello,

I want to use GO enrichment (http://geneontology.org/) for the first time with my E.coli proteomic hits. My IDs are from Uniprot and I have 211 of these. After submitting my 211 IDs, only 4 are mapped.

The 4 IDs that mapped successfully are:

```
P39356
P0A8A0
P0AAB8
P0AAI3
```

Here is a sample list of IDs that didn't match:

```
P0A954
A0A061KCV1
P0AGK6
A0A066SYL5
```

I have also tried to convert the Uniprot IDs to Uniprot Gene names and then searching it on GO, which has given me an improved number of mapped IDs to 129 out of 211.

A sample of hits which was mapped successfully includes:

```
ampC
accA
yeaG
atpF
```

Here are a sample of Gene names which was not mapped successfully:

```
bla
groL1
hlyA
AC789_1c11180
BY96_12590
```

If someone can shed some light on why it is not mapping successfully that will be great. Thanks.

The list of IDs I want submitting are located here https://gist.github.com/anonymous/01877dec923f6ce6f777

Answer: If the genome of your E.coli strain is publicly available you can get the GO annotation from QuickGO. Basically download the annotation of your stain and then extract the genes in your list.

Otherwise you could try to annotate your genes with Blast2GO or run them with Interproscan (there is also a stand alone version for the last one).

For making the enrichment there are many nice tools. I personally like the one in Bioconductor R, the most (e.g. topGO).

You can get an idea checking previous posts:

- https://www.biostars.org/p/153475/
- https://www.biostars.org/p/100223/
- https://www.biostars.org/p/128809/#128811
Question: I am trying to work out what tools use the secondary alignment feature of the sam format. My understanding is that there will always be exactly 1 primary alignment for each read, and secondary alignments could arise because the aligner cannot decide the exact best alignment.

Does BWA / other aligners emit secondary alignments?

Which downstream tools use secondary alignments? Haplotypecaller?

Answer: For RNAseq quantification, RSEM, eXpress and Salmon (when quantifying reads from BAMs) include secondary alignments in their quantification. Probably all mappers output secondary alignments, but implementation details may vary, and often there are flags to alter this behaviour. For example, `bwa mem -M` will mark shorter split hits as secondary instead of suplementary, or `STAR --outSAMprimaryFlag AllBestScore` will mark all multiple alignments with the same best score as primary.

I am not aware of other tools that use secondary alignments, but this doesn't mean they don't exist.
Question: In attempting not to reinvent the wheel, are there any lightweight tools that take a single BED12 element and quickly render its intron/UTR/CDS structure and directionality?

I'm not looking for a heavyweight toolkit or framework that renders a static genome browser shot of multiple BED12 elements, but something that takes in an element and simply shows a rough sketch of its structure. Thanks for any tips!

Answer: a very basic html form converting a bed12 entry to SVG :

https://gist.github.com/lindenb/3dd99f6429b7dc2092e0d1c94d74ed65
Question: Hello everyone,

I was using Porechop (as everyone I guess) in order to remove adaptaters from ONT reads. This tools is now unsuported (https://www.biostars.org/p/282767/#282775 , https://github.com/rrwick/Porechop), and with rapid changes in ONT softwares, I'm really afraid it will soon be impossible to use this tools anymore, and I'm looking for alternativ tools.

I found the new tools of Ryan Wick Deepbiner ( https://github.com/rrwick/Deepbinner ), which is presented by him as the sucessor of Porechop, but it is unclear if the tools only demultiplexes reads or also removes the adaptaters from them. For example, I feel like it's impossible to use on non barcoded runs.

Does anyone knows more about that subject ? I've been looking for some ideas on Nanopore Community, but so far my question remains.

Cheers,

Rox

Answer: Porechop may be no longer maintained, but it still works. I don't think there is any reason to quickly start looking for alternatives. Based on what I read on the nanopore community forum the adapter removal step would be added to Guppy, the basecaller which is to replace albacore in the "near future". So they (finally) start supporting this themselves.