cannin/biostars_qa huggingface dataset:
Question: Hi everyone:
I am interested in quantifying change in repetitive elements ( LTR here) transcription after treatment and I come up with following ideas:

1. Directly map RNA-seq data to genome with hisat2 and quantify with repetitive element annotation from Repeatmasker, followed by collecting elements from the same class to compare them. But I am not sure about how to set up maximum allowed multiple alignment value (For most RNA-seq it requires to be uniquely mapped but the value would be much higher since repetitive elements happens lots of times).
2. I got consensus repetitive element sequence fastq from Repbase, is it possible to view these repeat elements as "transcriptome" and use salmon (or similar transcriptome based tools) to map reads on it?

I am not familiar with this area and I would appreciate any suggestions . Thanks for help!

Update: Since I am only interested in LTR, I have modified the question. It looks possible to extract uniquely mapped reads and combine with Repeatmasker annotation. Direct quantification looks like will fail since repetitive elements are abundant in mRNA.

Answer: You can use STAR and then put though [TEtranscript](https://github.com/mhammell-laboratory/tetoolkit). You can allow multiple entries with STAR and it generally produces better alignments than hisat2 (in my experience at least). Our group that works on repeat elements uses this method.

While you can use the consensus repeat sequence, you end up biasing things for how close the expressed repeats are to the consensus. Consensus sequences are mostly useful for showing a profile over a single instance where you can label structure easily.
Question: 2 questions:

1. When use RNA-Seq data to do Differentially Expressed Genes (DEG) analysis, should the sample (/replicate) numbers of two groups must be the same? For example, if I have 8 samples from control, and 5 samples from treatment group, is it OK to use DESeq to do DEG analysis?

2. I'm using HISAT2 and featureCouts, after that, got the counts files, before putting them into DESeq, should I do normalisation firstly or can I use them directly?

Answer: I am expecting you are using DESeq2 and not DESeq.

- As far as numbers in each group is concerned, it is pretty fine to
perform DE analysis. The ideal scenario you get equal samples that
are paired and you need to use that feature while performing DE
analysis with any standard tool like DESeq2, edgeR or Limma.
- DESeq2 can still perform DE analysis with just 2 samples in one group
and 3 in the other. That's the lowest limit, going lower than that
the results are usually not trusted worthy.

- I have read edgeR can do even with lesser samples in the group but I
do not trust such analysis tbh. Your number of samples per group is
pretty good to perform the tests.
- About the normalization. The DE tools I mentioned and also you put in
query work on count data. So there is no point of putting normalized
data in them. They will perform normalization in the subsequent
steps. Just prepare your count table well and follow the DESeq2
tutorial and you are good to go.
- I will advice to follow the tutorial pretty well before performing
any DE analysis. It is always good to understand how the data behaves,
not only a QC ploy but also a good practice in exploratory data analysis.
Gives an understanding why you need to perform DE analysis and which samples
should be included in it. There might be a scenario where you might
have to move 1 or 2 samples from either of the group if they behave
as outliers, owing to either batch or sequencing errors, even if you
take care of them using batch correction methods. So a complete
workflow is advised and also enables you to make a pipleline for discovery which you might be using more often in your lab setting. I hope this was informative for your query.
Question: Dear all,

I think I’ve collected most of the available posts

from Biostars about RNA-secondary structure prediction software:

(T-boxes, riboswitches and attenuators search).

Most of them are 5-6 years old, but some are recent.

I hope someone else has made a better collection,

simply I didn’t find it. Please, help me!

This one, https://www.biostars.org/p/78883/#78888

is very nice, but 5.3 years old.

Search in NCBI, for example, gives more than 2000

articles, but it’s impossible to look through all of them.

Thank you!!

Natasha

RNA_secondary structure:

https://www.biostars.org/p/126368/

https://www.biostars.org/p/119924/

https://www.biostars.org/p/18141/

https://www.biostars.org/p/79426/

pictures

https://www.biostars.org/p/15569/

https://www.biostars.org/p/316248/

hairpin, attenuators

https://www.biostars.org/p/189815/

https://www.biostars.org/p/265768/#286744

https://go.mail.ru/search?fm=1&q=search%20for%20attenuator%20in%20RNA%20secondary%20structure&sbmt=1542733989299&src=go

search for attenuator in RNA secondary structure

https://www.biostars.org/p/274611/#274650

https://www.tbi.univie.ac.at/papers/Abstracts/martin_diss.pdf

Answer: RNA_secondary structure - (2017) http://www.bioinf.uni-freiburg.de/Publications/Fallmann_Will_Engelhardt-Recen_advan_RNA-2017.pdf
Question: Hello All,

I am trying to understand the sam output from the alignment. following is example:

[query_name]	16	[Chr_name]	[position] 28	[cigar]	*	0	0	[long_seq]	*	NM:i:11017	ms:i:912	AS:i:882	nn:i:0	ts:A:+	tp:A:P	cm:i:170	s1:i:835	s2:i:0

from [This link][1] I understand the what the first 11 field means. However, I am unable to understand the fields after 11th fields. Can someone help me understand them?

Thank you in advance.

[1]: https://samtools.github.io/hts-specs/SAMv1.pdf

Answer: In the link you gave, look at section 1.5 "The alignment section: optional fields." It says

> All optional fields follow the TAG:TYPE:VALUE format where TAG is a
> two-character string that matches /[A-Za-z][A-Za-z0-9]/

And it explains what those tags, types, and values are.
Question: Hi,
please, how can I visualize nucelotides for every position. I have a percentage of bases but now I don't know how to create a plot.

fastq <- readDNAStringSet("https://d28rh4a8wq0iu5.cloudfront.net/ads1/data/ERR037900_1.first1000.fastq","fastq")
fastq
freq <- alphabetFrequency(fastq, as.prob = T,baseOnly=T)

I would like something like this (sorry for quality):

![][1]

[1]: https://image.ibb.co/cQPULG/base.png

Answer: library(Biostrings)
fastq <- readDNAStringSet("https://d28rh4a8wq0iu5.cloudfront.net/ads1/data/ERR037900_1.first1000.fastq","fastq")
freq <- alphabetFrequency(fastq, as.prob = T,baseOnly=T)
matplot(freq,type='l')
legend(legend = colnames(freq),x=700,y=0.7,lty=1:5,col=1:5)

Note answer about consensusMatrix. That code is the best answer.
Question: I want to compare both k means and k harmonic means algorithm in terms of cluster quality. Can anyone suggest a good measure for comparison. Thanks in advance.

Answer: Most cluster metrics will be some form of inter. vs intra. cluster distances. The gist of these metrics is to see how how close are members within a cluster and how different are members among clusters. There really isn't a single good metric.

It depends on the input data and what you are looking for in clustering. The best way to compare cluster algorithms is to use real data for benchmarking.
Question: Hello everyone!
I would appreciate some help with this data frame 'df':
# df
Marker Sample allele
FlaHDF11	CP26	h
FlaHDF11	CP26	a
FlaHDF12	CP26	e
FlaHDF12	CP26	f
FlaHDF11	CP27	g
FlaHDF11	CP27	h
FlaHDF12	CP27	t
FlaHDF12	CP27	z
I would like something tlike this:
FlaHDF11	FlaHDF11	FlaHDF12	FlaHDF12
CP26	h	a	e	f
CP27	g	h	t	z

#This is the code Im using and I suspect the problem is in the line (6). Im always expecting to have just two rows for the same sample, but the code is optimized to add extra columns in case some samples have more than 2 alleles.

a = read.table('df', stringsAsFactors=F,header=T)
a = a[-1,]
a1 = a[!duplicated(a[,1:2]),]
rownames(a1) = paste(a1$Marker,a1$Sample)
t=1;
while(sum(duplicated(a[,1:2])) >0) #### I think the problem is here
{
a1 = cbind(a1,rep(NA, dim(a1)[[1]]))
a = a[duplicated(a[,1:2]),]
a2 = a[!duplicated(a[,1:2]),]
a1[paste(a2$Marker,a2$Sample),t+3] = a2$Size
t = t+1
}
colnames(a1) = c(colnames(a), paste('Size',1:(t-1),sep='.') )
a2 = a1[!is.na(a1$Size.2),]
m = matrix(NA,length(unique(a1$Sample)), length(unique(a1$Marker))*t)
rownames(m) = unique(a1$Sample)
colnames(m) = paste(rep(unique(a1$Marker), each=t), rep(1:4,length(unique(a1$Marker))), sep='.')
for (i in rownames(m))
{
tt = a1[a1$Sample %in% i,]
m[ i, paste(rep(tt$Marker,each=4), rep(1:4, dim(tt)[[1]]), sep='.')] = as.vector(t(tt[,3:6]))
}
m = m[,colSums(is.na(m)) < dim(m)[[1]]]

#### Error in `[<-.data.frame`(`*tmp*`, paste(a2$Marker, a2$Sample), t + 3, :
replacement has length zero

Kind regards!
Roberto

Answer: If we go ahead and "codify the assumption that there are two alleles" by renaming the marker column, then you can get the desired output similar to this

library(reshape2)
x=read.table(text='Marker Sample allele
FlaHDF11 CP26 h
FlaHDF11 CP26 a
FlaHDF12 CP26 e
FlaHDF12 CP26 f
FlaHDF11 CP27 g
FlaHDF11 CP27 h
FlaHDF12 CP27 t
FlaHDF12 CP27 z',header=T,stringsAsFactors=F)

x[seq(1,nrow(x),by=2),]$Marker=paste0(x[seq(1,nrow(x),by=2),]$Marker,'a')
x[seq(2,nrow(x),by=2),]$Marker=paste0(x[seq(2,nrow(x),by=2),]$Marker,'b')

acast(x,Sample~Marker)

#Output
#FlaHDF11a FlaHDF11b FlaHDF12a FlaHDF12b
#CP26 "h" "a" "e" "f"
#CP27 "g" "h" "t" "z"
Question: Hello！I'm using a package from github. I've installed it. But when i follow the guideline. I met a error about opening file.

require(rmotifx)
rmotifx

fg.path = system.file("extdata", "fg-data-ck2.txt", package="motifx")

bg.path = system.file("extdata", "bg-data-serine.txt", package="motifx")

fg.seqs = readLines(fg.path)

Warning message:
In file(con, "r") :
file("") only supports open = "w+" and open = "w+b": using the former

bg.seqs = readLines(bg.path)

Warning message:
In file(con, "r") :
file("") only supports open = "w+" and open = "w+b": using the former

head(fg.seqs)

character(0)

head(bg.seqs)

character(0)

Those two file are peptide sequence file. So i need to readLines them. I'm wondering how can i fix it?
Thank you!

here's the original tutorial: https://github.com/omarwagih/rmotifx

Answer: I believe this is an issue with how you are calling your data, rather than something with the `readLines()` function. After trying the example on my console I got the same error. The `system.files()` calls don't work so the objects `bg.path` and `fg.path` are empty, and when `readLines()` tries to search for the .txt files it can't find them. So, you just have to tell R exactly where to find the files you want it to run. You can do this by just including the entire file path in the `readLines()` call:

fg.seqs = readLines(file("C:/Users/Gregory/Documents/R/win-library/3.3/rmotifx/extdata/fg-data-ck2.txt"))
bg.seqs = readLines(file("C:/Users/Gregory/Documents/R/win-library/3.3/rmotifx/extdata/bg-data-serine.txt"))
head(fg.seqs)
###[1] "KRAGGEESQFEMDI_" "SSYGISETTLEEIFL" "SYGISETTLEEIFLK"
###[4] "FLKVAEESGVDAETS" "WSLNKEDTSEQVVPV" "KKLSVPTSDEEDEVP"

*EDIT: You'll need to find the files on your computer and change the path.
Question: I use the conda update the R language and then I require packages but there is a error.

> unable to load shared object
> '/home/slee_pkuhpc/lustre1/test/software/miniconda2/lib/R/modules//internet.so':
> libssl.so.1.0.0: cannot open shared object file: No such file or
> directory

I Googled the question but I did not find a fix. I do not have the root authority and do not change something.

So any advice is helpful.

Thanks.

Answer: conda install openssl==1.0.2p

What you're observing happens when you use `conda install` inside an environment that has packages with pinned dependencies for which an updated version of the dependency is available. Since `conda install` seemingly ignores pinnings, it gleefully updates everything in the environment thereby randomly breaking things. The best practices are one of the following:

* Use `conda create` to make a new environment for software packages/combinations
* Use `--no-update-deps` so conda doesn't try to update everything if you MUST use `conda install`
Question: Hi all,
I have an RNAseq dataset I am analyzing of gene expression from patients with a variety of illnesses. The design matrix is generally as follows:

> patient <- 1:12
> illness <- c(rep("control",3),rep("viral",3),rep("bacterial",3),rep("bactrial_and_viral",3))
> design <- model.matrix(~0+illness)
> rownames(design) <- patient

I want to find genes which are DE only in the bacterial group (relative to control), genes which are DE only in the viral group (relative to control), and genes which are DE only in the combined bacterial+viral group (relative to control). I am having trouble deciding which contrasts are appropriate to do this.

I know that if i just wanted to find genes which change in bacterial but not viral I could do:

> makeContrasts(bacterial_only = (bacterial-control)-(viral-control), levels = design)

But I think this would not exclude genes which also change in the bacterial_and_viral condition.

I cannot find the appropriate way to do the analogous comparison with more than 1 other experimental group. Could I do the following?

> makeContrasts(bacterial_only = (bacterial-control) - ( (viral-control) + (bacterial_and_viral-control) )/2, levels = design)

I'd really appreciate any insight on how to find genes which change in **only** one experimental group relative to control. I'm also worried I might be overthinking this--could this be accomplished by simply:

> makeContrasts (bacterial_only = bacterial-control, levels = design)

Does this find genes only DE in bacterial and not in other experimental groups?

Thanks for any insight you all can provide!

Brian

Answer: > I know that if i just wanted to find genes which change in bacterial but not viral I could do:
`makeContrasts(bacterial_only = (bacterial-control)-(viral-control), levels = design)`

That's not what the contrast does. The interaction contrast finds all genes for which the bacterial effect is significantly different to the viral effect. There's no assumption that a viral effect is absent. There might be strong effects in both bacteria and virus but in different directions. Or the effect could be just in bacteria. Or the effect could be just in virus. Or the effects in bacteria and virus could be in the same direction, but much bigger for one than the other.

You actually have only one control in this experiment, rather than separate controls for bacteria and virus. So the contrast
```
(bacterial-control)-(viral-control)
```
is exactly the same as
```
bacterial-viral
```
because the control cancels out. As you can see, the contrast just gives a differential expression test between bacterial and viral without regard to the control.
Question: I am looking for a software that can take SNPs information as input (e.g. a 23andme genomic data file) and haplogroup as output, and ideally with an explanation of the haplogroup.

If such a software does not exist, references to algorithms + databases would be also highly appreciated!

Answer: for mtDNA you can use HaploGrep 2 http://haplogrep.uibk.ac.at/ or Phy-Mer https://github.com/MEEIBioinformaticsCenter/phy-mer

for Y chromosome yHaplo https://github.com/23andMe/yhaplo
Question: Hi all,

I guess most people here are using open source/command line tools for NGS analysis. Nonetheless I am wondering if anyone has experience with any of the major commercial analysis suites and what your opinions are on them.

Undoubtedly there is a loss of flexibility and transparency when we move to an analysis suite but they have their obvious advantages also.

I am not so much concerned with proprietary algorithms etc so much as I am thinking of usability, stability, flexibility, etc...

Answer: I have been very impressed by the CLC Genomics Workbench and I have been recommending it to most of the life scientists regardless of their computational background.

It is particularly handy to study small genomes (10-100Megabases) that can be processed easily on a more powerful desktop. For many tasks I find it amazingly potent way to save time - rule out certain hypothesis, explore potential alternatives.

For example right now I have a dataset of 9 million 70bp reads on a bacteria of 4 Megabases. I can assemble that with three four clicks, or align to one reference or another, or blast the largest contigs - each task takes just a a few clicks and probably about five minutes of my time to set all that up - then I can do something else while it runs. If something good comes up I can move onto a more thorough process. Most impressively it does both the assembly and mapping faster and using less resources than most alternatives.

To use it I don't need to remember how to run Velvet or some other assembler, how to index the genome for mapping, pick the right parameters, etc. Compared to the overall cost of sequencing or that of my time the a license for CLC is very cheap.

If everyone could just use a tool like CLC and only asked for questions that cannot be solved by it - most people's research would be a lot farther ahead it would allow people to spend more of their brain power on complex tasks rather than simple things.
Question: Hello,

Does anyone know where I can find a bed file including repeat mask regions for the human genome?

Thanks!

Answer: Use [UCSC Table browser][1]. "Repeats" group and "Repeatmasker" track. Select output format as BED.

[1]: http://genome.ucsc.edu/cgi-bin/hgTables
Question: Dear community,

To 'flatten' a GTF is required for some programs. However, how does one flatten a GTF? There is not too much material online about it.

I had found this python script --written for Python 2.x, I believe-- for DEXSeq: https://github.com/olgabot/rna-seq-diff-exprn/blob/master/scripts/external/dexseq_prepare_annotation.py

However, is there a more rudimentary way?

Mit freundlichen Grüßen,

Kevin

Answer: Rsubread::flattenGTF()
Question: Hi all.

I am studying bioinformatics for my research and I have some questions regarding RNA seq vs. Microarray for gene expression level.

When considering their pros and cons, I think that RNA-seq has more good advantages over Microarray. For example, RNA is available for all species and also can check novel transcribed regions because RNA-seq platform doesn't require pre-defined procedures for checking gene expression level.

And also RNA has good advantages which are alternative splicing and to identify fusion genes. However, I didn't get the reason why RNA-seq has such two good points (alternative splicing, fusion genes) over Microarray techniques.

Could you explain the reason for me?

Answer: The two things you aren't sure about are due to microarrays only being able to detect exactly what they were designed for. In the case of fusion genes, to detect that on a microarray you'd have to have a probe spanning the fusion point. Similarly, for alternative splicing, you can only detect what you've designed probes for. That's one of the biggest gains with RNAseq, you can find things that you didn't explicitly want to look at beforehand.

I should note that the only real downside to RNAseq is that signals from genes/transcripts/etc. are competitive. So if you ever want to deconvolve signals originating from multiple sources (e.g., you have heterogenous samples and are interested in differential expression due to something within each of the sources rather than simply between treatment groups) then this is simpler with microarrays (i.e., things like independent component analysis are simpler).
Question: I was just wondering how does sliding window in trimmomatic work?

The definition is scanning from the 5' end of the read, and removes the 3' end of the read when the average quality of a group of bases drops below a specified threshold. For example, if we have a sequence `ATCGATCGATCG` and we set `SLIDINGWINDOW: 4:15`.

It begins with the first 4 in a window, `ATCG`, but if the score is below 15, which base it will trim? Is that the last base in this window? What is the next start position of the window? 2 or 5? thanks.

Answer: The SLIDINGWINDOW trimmer will cut the leftmost position in the window where the average quality drops below the threshold and remove the rest of the read. However if there is low quality in the very beginning of the read then it will fail the minimum length tests and be removed completely - the remaining 3-prime end (even if it is good quality will not be printed)

Consider the following test file `t.fq`:

@1\1
AATGATCGTAGCGATGCAAGCTAGCCCGATGCCCGATCGCATCG
+
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeEFCB
@2\1
AATGATCGTAGCGATGCAAGCTAGCCCGATGCCCGATCGCATCG
+
EFCBeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee

and processing with:

$ java -jar trimmomatic.jar SE -phred64 t.fq tt.fq SLIDINGWINDOW:4:15

TrimmomaticSE: Started with arguments: -phred64 t.fq tt.fq SLIDINGWINDOW:4:15
Automatically using 16 threads
Input Reads: 2 Surviving: 1 (50.00%) Dropped: 1 (50.00%)
TrimmomaticSE: Completed successfully

The output file looks like the following:

@1\1
AATGATCGTAGCGATGCAAGCTAGCCCGATGCCCGATCGC
+
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee

As you can see the read with the poor quality at the beginning has been removed completely
Question: After aligning my raw RNA-Seq reads with Salmon, I imported the quant files into R and created a count matrix with annotations using tximeta (creating a SummarizedExperiment object).

Afterwards I built a DGEList object in order to pass it onto edgeR, using the tximeta function makeDGEList. This creates an offset matrix, which I believe normalizes with respect to average transcript length.

Now, when I tried normalize the data in edgeR using calcNormFactors, it gives me this warning:

LL_normalized <- calcNormFactors(LL_genefiltered, method = "TMM")

>Warning message:
In calcNormFactors.DGEList(LL_genefiltered, method = "TMM") :
object contains offsets, which take precedence over library
sizes and norm factors (and which will not be recomputed).

It appears the offset matrix from tximeta and the normalization attempted by edgeR conflict with one another, and they have different norm factor values.

I would like to normalize by sequencing depth, RNA composition (effective library size), and gene length.

So my question is, does the offset matrix created by tximeta's makeDGEList function account for all these factors (I don't think it does, only for gene length), or do I need to pass on the DGEList object without an offset matrix and let edgeR take care of the normalization(if so, how?), or is there a way to make them work together?

Of note, the tximport vignette warns not to manually pass the original gene-level counts to downstream methods without an offset. Perhaps this was assuming that no other normalization would be done after tximport/tximeta, and can be disregarded?

Answer: The offset matrix includes both correction for average transcript length **and** library size / composition.

See here the code: https://github.com/mikelove/tximeta/blob/master/R/dgelist.R#L12
It includes the `calcNormFactors` already.

Here is a more commented version what the steps do at tximport vignette: https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html#edgeR

That means you don't have to do anything towards normalization, the DGEList from `tximeta` is ready to go. The next logical step for edgeR would be `filterByExpr` and then the DE testing routine as described in the edgeR manual.
Question: What methylation values are usually supplied within ncbi geo 450k series matrices?

I am going to automatically download a lot of GSE (~80) matrices and to use them in my analysis. Before this I want to double check here if I understand correctly what data is within those matrices.

I know that they contain sample, contact information, etc. But what about methylation data? Are those values already **normalised, batch corrected**? Or does this vary between each study?

Answer: The methylation values are included in the GSE, yes. The values vary between studies. In general, they will be normalized, but the normalization methods will vary from GSE to GSE.