cannin/biostars_qa huggingface dataset:
Question: Hi there,

I know that this is probably a common and newbie question but I can't find the solution. I have a strand specific RNA-seq data, specifically first-strand. Which is the correct way to specify this in trinity, RF or FR? I've read the manual but is not clear for me...

Thanks in advance.

Answer: The Trinity documentation says:

If you have strand-specific data, specify the library type. There are four library types:

- Paired reads:
- **RF**: first read (/1) of fragment pair is sequenced as anti-sense (reverse(**R**)), and second read (/2) is in the sense strand (forward(**F**)); typical of the dUTP/UDG sequencing method.
- **FR**: first read (/1) of fragment pair is sequenced as sense (forward), and second read (/2) is in the antisense strand (reverse)
- Unpaired (single) reads:
- **F**: the single read is in the sense (forward) orientation
- **R**: the single read is in the antisense (reverse) orientation

The [TopHat manual](https://ccb.jhu.edu/software/tophat/manual.shtml) says:

Library Type Examples Description
fr-unstranded Standard Illumina Reads from the left-most end of the fragment (in transcript coordinates) map to the transcript strand, and the right-most end maps to the opposite strand.
fr-firststrand dUTP, NSR, NNSR Same as above except we enforce the rule that the right-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during first strand synthesis is sequenced.
fr-secondstrand Ligation, Standard SOLiD Same as above except we enforce the rule that the left-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during second strand synthesis is sequenced.

To complete that, in Trinity `fr-firststrand` corresponds to `RF` and `fr-secondstrand` corresponds to `FR`.

You can have a look here: http://rnaseq.uoregon.edu. They well explain the difference between first-strand and second-strand synthesis.

You can also have look to this publication: Sequencing technologies - the next generation. Nature reviews. Genetics. 2010. doi:10.1038/nrg2626

Briefly, if you know the technology used for the sequencing you should be able to guess.
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: R packages such as APE and Phangorn and so on allow you to compute evolutionary distance measures. Both packages and other related, phylogenetic based packages are written - at least in part - by the same author: Paradis.

The function dist.dna, for example, accepts as input an aligned set of DNA sequences, and has an optional parameter "model", which allows you to choose from many different evolutionary models, including all the ones you are looking to use here. http://www.inside-r.org/packages/cran/ape/docs/dist.dna. If you'd like to see how that function implements various models you could take a look at the code - although some of it is likely written in C as well as R.

If you read a book like the Phylogenetic Handbook [here][1], where they discuss the mathematical derivations of distance measures - in fact there is a chapter dedicated to it - you will learn the classic models are in fact specific cases of the GTR model.

[1]: http://www.cambridge.org/gb/academic/subjects/life-sciences/genomics-bioinformatics-and-systems-biology/phylogenetic-handbook-practical-approach-phylogenetic-analysis-and-hypothesis-testing-2nd-edition
Question: Hi,

I am working on the expression datasets that has multiple-variables. I am using the `FactoMineR` and `pca3d` libraries for this purpose and was able to distinguish two factor levels belonging to a one variable column. However, I am not able to perform the PCA on the two variable columns that has different factor levels. Please let me know how can I do the PCA on multi-variable.

Below is the code I have ran with the one variable column;

Neg_Dct
## Column [1-4 and 269-272] of the data set contains variables/categorical data##
df = Neg_Dct[,-c(1:4, 269:272)]

library(FactoMineR)
nb_1 = estim_ncpPCA(df,ncp.max=5)
res.comp_1 = imputePCA(df,ncp=2)
res.pca_1 = PCA(res.comp_1$completeObs)

library(pca3d)
Node <- res.comp_1$completeObs
pca <- prcomp(Node, scale.=TRUE)
gr <- factor(Neg_Dct[,272])
summary(gr)
B_1 B_2
31 132

#2D plot##
pca3d(pca, group=gr)

#2D plot##
pca2d(pca, group=gr)

Thank you,

Toufiq

Answer: Use PCAtools (I am the developer), where you can have as many levels as you want, and represent them as both different shapes and / or colours: 4.2.3 Change shape based on tumour grade, remove connectors, and add titles
Question: Hi,
I am trying to convert a Canine gene annotation (GTF) file downloaded from Ensembl to BED file using the gtf2bed tool within the BEDOPS application. Using this command gives an error:

$ gtf2bed < Canis_familiaris.CanFam3.1.85_noheader.gtf > Canis_familiaris.CanFam3.1.85_noheader.bed

Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)

I checked the first few lines of the GTF file and it seems to match up with the required format:

$ head Canis_familiaris.CanFam3.1.85_noheader.gtf

X ensembl gene 1575 5716 . + . gene_id "ENSCAFG00000010935"; gene_version "3"; gene_source "ensembl"; gene_biotype "protein_coding";

X ensembl transcript 1575 5716 . + . gene_id "ENSCAFG00000010935"; gene_version "3"; transcript_id "ENSCAFT00000017396"; transcript_version "3"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding";

...

I looked at the source code on github for this tool and can see that is check for gene or transcript id and if not present gives this error. But the gene_id is present here in the first line, so not sure how it is reaching the error condition.

I would appreciate any help with troubleshooting this error.

Thank you,
- Pankaj

Answer: **Note: This should no longer be an error with *gtf2bed* v2.4.40 and on.**

----------

I added more stringent GTF format validation to BEDOPS v2.4.20.

The error suggests that the first line is missing the `transcript_id` field. It has a `gene_id` field, as you note, but no `transcript_id` field. The GTF 2.2 specification indicates that this field is mandatory, though its value can be an empty string.

There are a couple solutions:

1. Use an older version of `gtf2bed` that doesn't apply this validation check (e.g., 2.4.19 or earlier)
2. Or, modify the GTF and add a placeholder field where none exists

I suggest the second solution. You could do the following:

$ awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' input.gtf | gtf2bed - > output.bed

This adds `transcript_id "";` to lines in the GTF that do not contain that field, and leaves other lines unchanged.

The GTF that comes out of this `awk` statement is more valid, enough to get through the conversion step, and so it can be piped to `gtf2bed` to get BED as output.
Question: Hi Biostars,

Is it legitimate to sum up raw read counts from technical replicates of RNAseq, and use these summed counts for DE analysis. Would appreciate detailed and justified answers.

Thanks

Answer: There is a theoretical justification for summing and not averaging. Read counts follow a Poisson distribution so averaging them results in data that is not Poisson, but summing is still Poisson. See Mike Love's explanation [here][1]

[1]: http://seqanswers.com/forums/showthread.php?t=60996#4
Question: Dear everyone

I'm confused that how to build Risk score formula for survival analysis. I have been seen many articles involved the Risk score formula building, but I have spent a day to solved it but failed. can someone taken a example to show how to build it?

[This is the article details][1]:

> Base on the expression level of five lncRNAs, we designed a risk-score formula for ccRCC patients’
survival prediction. The risk score formula is as following:

Risk score= (1.43 × expression level of AC069513.4) +
(0.81 × expression level of AC003092.1) +
(1.64 × expression level of RP11-507K2.3) +
(-6.56 × expression level of CTC-205M6.2) +
(-1.72 × expression level of U91328.21)

but I don't know how to calculate the Risk score formular

Thank you very much!

Alex

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5601685/

Answer: As to which zx8754 alluded, it is simply about multiplying the expression values by a weight, and then summing these up.

The weights that they use are the beta estimates (coefficients) from the fitted Cox proportional hazards model. The exponential function applied to these coefficients gives the odds ratios.

This, and variations of it, represent a very common way to calculate 'risk scores'.

Hopefully that helps.

Kevin
Question: Simple question - I need to create a GTF file to use in HTSeq-count that contains gene regions plus 3kb upstream. (Background: doing a MeDIP-seq experiment, want to look for differential methylation in genic and 3kb promoter regions using count based method like edgeR/DESeq).

I was planning on making one myself from the UCSC hg19 refFlat table. The refFlat table has gene coordinates, but I need to extend this 3kb upstream to capture promoter regions.

Column 3 contains the strand (+/-) and columns 4-5 contain the transcription start (txStart) and end (txEnd) positions.

If I want to capture 3kb upstream of the TSS, I was planning on adding 3000 to txStart, but is only for genes on the + strand, correct? If I want 3kb upstream of the TSS for genes on the - strand, should I add this 3000 to txEnd?

E.g. the DENND1B gene is on the - strand at chr1:197,473,879-197,744,623. However, looking at it in the browser, it's transcribed "right to left", so I presume I would want to add 3000 to the txEnd number, 197,744,623, even though this is really where transcription starts.

Am I thinking about this correctly?

Answer: > "If I want 3kb upstream of the TSS for genes on the - strand, should I add this 3000 to txEnd?"

yes.

the position are always ordered from the 5' of the chromosome to the 3' and named xxStart and xxEnd (with xxStart<=xxEnd) whatever the value of "strand".

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e 'select count(*) from refFlat where txStart>txEnd or cdsStart>cdsEnd'
+----------+
| count(*) |
+----------+
| 0 |
+----------+

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e 'select txStart,txEnd from refFlat where strand="-" limit 4'
+-----------+-----------+
| txStart | txEnd |
+-----------+-----------+
| 62929370 | 62937380 |
| 34610 | 36081 |
| 153270337 | 153283194 |
| 113235158 | 114449242 |
+-----------+-----------+

So you'll only have to :

newStart=(start<3000?0:start-3000);
newEnd=end+3000;

**EDIT:** I was too fast. My pseudo-code extends the sequence on both side. If you want to extend **upstream** sequence :

newStart=(strand=='+'?start<3000?0:start-3000:start);
newEnd=(strand=='+'?end:end+3000);

to fetch the (positive strand) sequence; See

[How to get the sequence of a genomic region from UCSC?][1]

and

[Extract sequence from the genome?][2]

[1]: /p/56/
[2]: /p/3137
Question: Hello all,
I don't know if this is the appropriate place to ask to this question, but I am currently analyzing a single cell rna-seq dataset with 2 conditions and 10 samples each (total of 20 samples). The samples are not technical replicates, they come from different humans each. The files, once downloaded, are in .H5 format. After scouring the web, I still cannot find a pipeline that works for me. Should I analyze each of the samples individually or merge them into one some way and proceed? I am trying to follow the Seurat pipeline line so any help on that note would be helpful. I have written this code:

h5_files<-list.files(pattern ="*.h5")
h5_read <- lapply(h5_files, Read10X_h5)
h5_seurat <- lapply(h5_read, CreateSeuratObject, min.cells=5, min.features=250, project="ccc")

here h5_seurat in the last line of code is a list with 20 Seurat objects (from the original .H5) files.
I want to follow the rest of the Seurat pipeline, but I don't know if should do it individually for each of the 20 samples or merge them somehow?

Thanks

Answer: Hello,

You can merge all the files into a Seurat object and continue the analysis with it. Here's an example similar to what you're looking for [scRNA-seq example][1].

Best,

Rafael

[1]: https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_01_qc.html
Question: Hi guys !

From a bam alignment file, I want to compute the ratio between the number of reads terminating and overlapping at all genomic positions (called *psi-ratio* in this figure).

![image description][1]

For the number of reads overlapping, it's easy (its basically the coverage, given by `samtools depth` for instance), but for the reads terminating I'm a bit lost. I guess this has to do with the POS field of the sam/bam alignment (1-based leftmost mapping POSition) but I can't go much further than that... Any help or resource will be appreciated !

Subsidiary question : what if the input data is paired-end and I want the ratio between the **FRAGMENTS** terminating and overlapping ?

For those interested in why I want to do that, the idea is very similar to the recently described [psi-seq][2] : I want to detect positions on RNAs that blocks the reaction of reverse transcription during the library preparation.

Thanks a lot !

[1]: http://i.imgur.com/WlQko3x.png
[2]: http://www.ncbi.nlm.nih.gov/pubmed/25219674

Answer: > for the reads terminating I'm a bit lost.

using **bioalcidae** : https://github.com/lindenb/jvarkit/wiki/BioAlcidae

```
$ java -jar dist/bioalcidae.jar input.bam \
-e 'while(iter.hasNext()) { var rec = iter.next(); if(rec.getReadUnmappedFlag()) continue; out.println(rec.getContig()+"\t"+rec.getAlignmentEnd()); }' | LC_ALL=C sort -k1,1 -k2,2n | LC_ALL=C uniq -c

1 rotavirus	52
1 rotavirus	56
1 rotavirus	64
1 rotavirus	66
1 rotavirus	67
2 rotavirus	68
1 rotavirus	69
6 rotavirus	70
8 rotavirus	71
5 rotavirus	72
5 rotavirus	73
7 rotavirus	74
4 rotavirus	75
6 rotavirus	76
7 rotavirus	77
5 rotavirus	78
8 rotavirus	79
3 rotavirus	80
10 rotavirus	81
4 rotavirus	82
6 rotavirus	83
7 rotavirus	84
4 rotavirus	85
3 rotavirus	86
10 rotavirus	87
5 rotavirus	88
6 rotavirus	89
4 rotavirus	90
9 rotavirus	91
7 rotavirus	92
4 rotavirus	93
1 rotavirus	94
5 rotavirus	95
5 rotavirus	96
6 rotavirus	97
7 rotavirus	98
3 rotavirus	99
4 rotavirus	100
(...)
```

you might use getUnclippedEnd instead of getAlignmentEnd if you want the unclipped alignments.
Question: Hi pals, I got some problems here.

I have an excel file with several genomic coordinates in the format "strand:start:end", and I also have access to my genome of interest in different formats (fna, gbk...). Now, what I want to do is to extract and obtain a file with all the sequences corresponding to my coordinates, I don't mind if the sequences have a header or a custom ID, I just want the sequences in the same order as I have them in my excel file. Is there any tool that can do this? It would be much appreciated since I have to extract more than 2000 sequences, and doing this manually is humanly impossible.

The next step is to perform a BLAST with another genome, also accessible in several formats. I'm using the blastn NCBI program and this script `-task blastn -query "query" -subject "subject" -oytfmt 6 >name`, how can I make the output file to only contain the most relevant matches? (E.g. only the matches that have an E value lower than 0.0001)

Thanks beforehand!

Answer: If you're using Windows, save your excel file as text with tab delimiters, and first column chromosome name, second start, third end (this is BED format. Then install Python3, then open `cmd.exe:`

```
C:\Python34\Scripts\pip.exe install pyfaidx
C:\Python34\Scripts\faidx.exe genome.fna --bed regions.bed > out.fasta
```
Question: I used makeblastdb to search many short fasta sequences in a known organism. I successfully completed this step and got the mappings. My question is

How can I get the start and stop coordinate?

for eg

Query 1 AATATAGGTGGTACCACGGAATATCCGTCCTATTTGTATATAGGATGGATAtttttattt 60
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 1765181 AATATAGGTGGTACCACGGAATATCCGTCCTATTTGTATATAGGATGGATATTTTTATTT 1765122

Query 61 ttttAGGAGGTATAGCAAATGG 82
||||||||||||||||||||||
Sbjct 1765121 TTTTAGGAGGTATAGCAAATGG 1765100

How to get 1765100 and 1765181 from a text file with many mappings like this. I would also like to count the number of mappings or query sequences?

Answer: blastn -query some.fa -db some.fa -outfmt 6 | awk '{print $9 "\t" $10}'
Question: I have a list of genomic ranges mapped to hg19 . my data is in matrix format lets call it `ranges` which has 600,000 rows and 4 clumns
here is few row of my data

head(ranges)
chr start end strand
[1,] "chr1" "10025" "10525" "."
[2,] "chr1" "13252" "13752" "."
[3,] "chr1" "16019" "16519" "."
[4,] "chr1" "96376" "96876" "."
[5,] "chr1" "115440" "115940" "."
[6,] "chr1" "235393" "235893" "."
Is there a function that gets sequences and calculates GC content for each row( each range)
I would prefer that output be in a vector format
I would really appreciate your help

Answer: A solution using GenomicRanges and BSgenome from Bioconductor. You will need to transform your data to a GRanges object for this.

library(GenomicRanges)
library(BSgenome.Hsapiens.UCSC.hg19)
library(BSgenome)

gr <- GenomicRanges::GRanges(seqnames = c("chr1", "chr2"),
ranges = IRanges(start = c(123456,123490),
end = c(500020, 600020)))

GetGC <- function(bsgenome, gr){

seqs <- BSgenome::getSeq(bsgenome, gr)
return(as.numeric(Biostrings::letterFrequency(x = seqs, letters = "GC", as.prob = TRUE)))

}

> gr
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand

[1] chr1 123456-500020 *
[2] chr2 123490-600020 *
-------
seqinfo: 2 sequences from an unspecified genome; no seqlengths

> GetGC(bsgenome = BSgenome.Hsapiens.UCSC.hg19, gr = gr)
[1] 0.2823072 0.4355687
Question: I am currently tinkering with an analysis bioconductor ppackage found here:

"An end to end workflow for differential gene expression using Affymetrix microarrays"
https://bioconductor.org/packages/devel/workflows/vignettes/maEndToEnd/inst/doc/MA-Workflow.html

Under section 4.5 you perform the task of running a general linearized model where sva surrogate variables are modeled as a function of your covariates. The following does this by pulling from an expression set object created in the first part of the pipeline:

glm.sv1 <-glm(pData(inpData_sv)[,"sv1"]~pData(inpData_sv)[,"Batch"]+pData(inpData_sv)[,"Sex"])

Unfortunately, I am getting an error telling me:

Error in pData(inpData_sv) : could not find function "pData"

I have tried getting the pData function by downloading ballgown and phylobase, the only two packages where I have seen pData(). Neither of these have seemed to provide me with the pData function required here. I'm pretty new to R, so could anyone recommend a substitute function for this R command?

Answer: From the linked manual:

> The *pData* function of the Biobase package...

Do you have Biobase package installed/loaded?
Question: How to calculate effect size between two samples for a gene set in case if I have effect size of the two samples?

for example

s.dist sample1 sample 2

0.4730	-0.05460 -0.47000

How s,dist has been calculated?
Thanks

Answer: > How s,dist has been calculated?

Since you have *euclidean* among your keywords, I suspect the answer is the euclidean distance between the two samples:

```
sqrt((-0.05460)^2 + (-0.47000)^2)
0.4731608
```

perhaps the discrepancy with your 0.4730 is due to some rounding error in reporting.
Question: I've a table of blast results with ~5hits/query protein:

```
Protein Class
ProtA 1
ProtA 1
ProtA 1
ProtA 0
ProtA 1
ProtB 1
ProtB 1
ProtB 0
ProtB 0
ProtB 1
```

I would like to convert this into a feature vector matrix like this:

```
Protein Class1 Class2 Class3 Class4 Class5
ProtA 1 1 1 0 1
ProtB 1 1 0 0 1
```

Can someone suggest me an efficient way to do this, since I've ~2300k hits in the file.

Answer: For any number of hits per protein:

```
> tmp
protein class
ProtA 1
ProtA 1
ProtA 1
ProtA 0
ProtA 1
ProtB 1
ProtB 1
ProtB 0
ProtB 0
ProtC 1
ProtD 1
ProtD 0

# install & load packages
library(reshape2)
library(plyr)
library(doMC)
registerDoMC(8) # assign 8 cores for ddply to use

# define function myfunc
myfunc <- function(x)
{
dim(x)
x$name = paste("Class",1:nrow(x))
dcast(x,protein~name,value.var = "class")
}

# call myfunc in ddply. run ddply in parallel mode to use >1 cores
res = ddply(.data = tmp,.variables = "protein",.fun = function(x) myfunc(x),.parallel=T)

# output
> res
protein Class 1 Class 2 Class 3 Class 4 Class 5
ProtA 1 1 1 0 1
ProtB 1 1 0 0 NA
ProtC 1 NA NA NA NA
ProtD 1 0 NA NA NA
```
Question: For each of several genomes, apart from the already available fasta sequence and associated GFF3 annotation files, I have also generated 5 additional GFF files for start-stop coordinates of 5 additional types of genomic features.

My goal is four-fold in the context of these 6 GFF files and 1 genomic DNA sequence.

1. **Explore** where two or more of these genomic features overlap / intersect / co-localize - I am doing this via text manipulation, using bedtools overlap or bedtools intersect. But I want to visualize 2 or more tracks , and just text-based calculation of intersection is not satisfying during exploration phase. I want to visualize it.

2. Generate **high-detail images** (with flexibility of color / shapes like that in IGB or gff2ps) for **small intervals**, looking at specific and most interesting cases.

3. Generate **overview images** for **entire** chromosomes or even a **genome**, for overlap across these 6 types of genomic features , without confusing or overwhelming the reader.

4. Finally, to request advice on which tool and/or which **statistical test** to perform for verifying whether the observed physical **co-localization** of any 2 of the 6 types of genomic features is random or **non-random**. And if latter is true, are there more sophisticated tests to examine physical distribution of genomic loci types, relative to one another? And are there tests that can examine more than 2 types of genomic features at a time?

The rather old thread at https://www.biostars.org/p/363/ discusses answers to questions 1 - 3 above, but I am curious to know if there are better / updates tools for my goals than the ones I mentioned above or at the link (bedtools or bedOps, IGB, GBrowse, GFF2PS). Thanks!

Answer: Hi,

If you can use R you should be able to create these plots with [karyoploteR](	http://bioconductor.org/packages/karyoploteR/). You would need to load the data into R (probably using [rtracklayer](http://bioconductor.org/packages/rtracklayer/)'s ' import' function) and then plot them using `kpPlotGenes`for genes and `kpPlotRegions` for everything else. You can find more information and various examples on how to use it at [karyoploteR titorial page](https://bernatgel.github.io/karyoploter_tutorial/).

As for point 4 you can use the Bioconductor package [regioneR](http://bioconductor.org/packages/release/bioc/html/regioneR.html). If you load the data into R you can use the function `overlapPermTest` to perform a permutation strategy to test if two sets of genomic regions overlap more (or less) than expected by chance.

Hope this helps

Bernat