cannin/biostars_qa huggingface dataset:
Question: Hi all,

I would like to get the upstream and downstream intergenic distance between two genes in R (taking into account that some genes may overlap). For example - in the following case I would expect to get a the upstream distance of 5 for gene A, a distance 0 or negative for downstream of A, and a distance of 10 for downstream of B/upstream of C.

-----|____gene A____|-------------------------|____gene_C___|--- + strand
---------------|______gene_B_______|---------------------------- - strand

I tried to use precede() and follow() in GRanges but these two functions completely exclude overlapping genes from the analysis. This sounds like a task everyone is doing and it should be part of a package but I can't find anything appropriate.

Thank you for any hints!

Answer: Here is some attempt. I create an example ranges (IRanges for simplicity). You can used this modified version of `plotRanges()` function (from the IRangesOverview vignette in the [IRanges][1] package) to visualize the ranges and how they overlap.

The problem using follow() and precede() is that they ignore overlapping ranges. Therefore here instead compute the sequence-wide comparison of distances. Maybe this needs to be modified to be more efficient and only do this in the *closest* sequence. However, with overlapping sequences it is not so clear to me what *closest* means.

library(IRanges)
foo <- IRanges(start = c(1, 7, 10, 5, 13), end = c(4, 12, 15, 8, 14))

plotRanges <- function(x,
xlim = x,
main = deparse(substitute(x)),
col = "black",
sep = 0.5,
...)

{
height <- 1
if (is(xlim, "Ranges"))
xlim <-
c(min(start(xlim)), max(end(xlim)))
bins <-
disjointBins(IRanges(start(x), end(x) + 1))
plot.new()
plot.window(xlim, c(0, max(bins) * (height + sep)))
ybottom <-
bins * (sep + height) - height
rect(start(x), ybottom, end(x), ybottom + height, col = col, ...)
text(c(start(x) + width(x)/2) - .5, y = ybottom + height /2)
title(main)
r <- range(x)
axis(1, at = seq(start(r), end(r)))
}

plotRanges(foo, col = rainbow(length(foo)))
foo
IRanges object with 5 ranges and 0 metadata columns:
start end width

[1] 1 4 4
[2] 7 12 6
[3] 10 15 6
[4] 5 8 4
[5] 13 14 2

I compute the sequence-wise comparison as this:

# comparison matrix.
combmatrix <- t(combn(length(foo), 2))
combmatrix
colnames(combmatrix) <- c("i", "j")
i j
[1,] 1 2
[2,] 1 3
[3,] 1 4
[4,] 1 5
[5,] 2 3
[6,] 2 4
[7,] 2 5
[8,] 3 4
[9,] 3 5

Then iterate over the comparison matrix and compute the distances:

hoo <- apply(combmatrix, 1, function(x) {
s1 <- foo[x[1], ]
s2 <- foo[x[2], ]
if (start(s1) < start(s2))
start(s2) - end(s1)
else
start(s1) - end(s2)
})

cbind(combmatrix, distance = hoo)
i j distance
[1,] 1 2 3
[2,] 1 3 6
[3,] 1 4 1
[4,] 1 5 9
[5,] 2 3 -2
[6,] 2 4 -1
[7,] 2 5 1
[8,] 3 4 2
[9,] 3 5 -2
[10,] 4 5 5

Note that for the comparison between sequences 3 and 5, it is not clear to me what the right answer would be in this case. I guess this is the reason distances to overlapping ranges are zero in IRanges (and GRanges).

[1]: http://bioconductor.org/packages/release/bioc/html/IRanges.html
Question: Hi all,

I am getting an error while running [the liftOver tool][1] locally on my Mac OSX

The code I execute is as follows:

~/Downloads/liftOver ./ICBP-summary-Nature.bed ~/Downloads/hg18ToHg19.over.chain.gz ./converted_icbp_coordinates.bed ./unlifted_icbp_coordinates.bed

The error reads as following:

```
Reading liftover chains
Mapping coordinates
Expecting integer field 2 line 2146483 of ./ICBP-summary-Nature.bed, got 6.5e+07
```

So I checked the data on this particular line and the data seems to be fine as follows by running sed as

sed -n '2146843p' ICBP-summary-Nature.bed

and I get the following output for the line

chr16 65819026 65819027 rs3852699

The consequence of this error is that the program breaks at this line number so my remaining snps (~300,000) in the subsequent lines do not get converted.

I thought I should post this query on the UCSC Genome Support Forum [here][2], but the registration seems to have closed?

Thank you for going through my query and I would be grateful for any tips or suggestions on how to resolve this?

**Update**:

At the request of komal, I am posting a few lines which appear before and after this line. So the following are some of the lines before 2146843(including 2146843) are:

```
chr16 65787923 65787924 rs3730406
chr16 65790185 65790186 rs11700
chr16 65791635 65791636 rs8058861
chr16 65791780 65791781 rs8059662
chr16 65794844 65794845 rs13336793
chr16 65798783 65798784 rs12051247
chr16 65798928 65798929 rs12051249
chr16 65802113 65802114 rs16957240
chr16 65811920 65811921 rs7193713
chr16 65817068 65817069 rs7196793
chr16 65817164 65817165 rs7196989
chr16 65819026 65819027 rs3852699
```

And the following are some of the lines afterwards:

```
chr16 65819026 65819027 rs3852699
chr16 65819348 65819349 rs6499116
chr16 65822861 65822862 rs6499118
chr16 65829359 65829360 rs3852700
chr16 65829878 65829879 rs16957265
chr16 65835167 65835168 rs6499119
chr16 65842341 65842342 rs6499121
chr16 65847342 65847343 rs9940665
chr16 65847656 65847657 rs9931407
chr16 65850656 65850657 rs8064216
chr16 65855664 65855665 rs8053031
```

[1]: http://genome.sph.umich.edu/wiki/LiftOver
[2]: http://redmine.soe.ucsc.edu/forum/index.php

Answer: Ok so I finally found out what was going wrong.

I misread line 2146483 as 2146843. Hence I could not find the error.

I split the file into two using the linux split command

split -l 1230663 ./Blood_pressure/ICBP-summary-Nature.bed ./Blood_pressure/split_blood_pressure

And that is when I realized my oversight as the 2nd file also gave the same error albeit in a different line number obviously which I fortunately read correctly!

Here is what the erroneous line in the file looked like:

chr16 6.5e+07 65000001 rs2289150

I changed it using sed as follows:

sed -i '' 's/6.5e+07/65000000/' split_blood_pressureab

Then I joined the two files through ´cat´and proceeded as usual with liftOver.

My sincere thanks to komal for taking the time to help me out.
Question: Dear all,

I have some questions regarding demultiplexing Illumina BCL files. First question is if Illumina´s Undetermined Fastq files contains barcodes?

Second question - If I would like to provide demultiplexing and keep my barcodes in fastq files (no in the header like CASAVA provide, but inside the read sequence) - does anybody have any experience how to do that?

Thak you for any sharing your experience.

Answer: MiSeq MSR does not put indexes in header but the index_rank (row). If the index used is not found in the sample sheet, it put a 0 (the last :0 in the header of the undetermined_reads.

With bcl2fastq you can ask to have the indexes written as a FASTQ with option `--create-fastq-for-index-reads

You can set MiseqReported to produce the fastq of the indexes modifying the MiSeqReporter.config with:

If you want to explore the indexes in the undetermined to see if you have mispelled an index in the sample sheet and you are loosing this sample in the demultiplexing, there is a file in the flowcell folder where are recorded all the index found an how many tiems they where found (you can acces this data also from the MSR) but I can not remember the path right now (this file is different in MiSeq and in NextSeq500).
Question: Hello

In Seurat there is a function to take the proportions of each cell identity so you can easily plot it with ggplots or something similar. However, most scRNA datasets I have seem (I mostly reanalyze data) have different sample sizes for each condition. So I'm sure just taking the proportions of cells might not be adequate. I believe you would need to normalize this. The first thing that comes to mind is dividing the number of cell identities by the number of conditions, but it still doesn't make much sense I guess, as sometimes the same conditions may have a high variation of cell identities too. [Here][1] the authors plot it by log2 of relative proportions, which I believe it is Z-score, but still it is a bit weird to me, as they have different numbers of samples in each status.

I couldn't find any Seurat vignette addressing this. Any solutions? Does my concern make sense?

[1]: https://immunology.sciencemag.org/content/immunology/5/49/eabd1554/F1.large.jpg?width=800&height=600&carousel=1

Answer: To compare cell proportions between conditions, I've found using a monte-carlo/permutation test to be the most sensible and robust way. The null hypothesis you want to test against is that the difference in cell proportions for each cluster between conditions is just a consequence of randomly sampling some number of cells for sequencing for each condition. To generate this null distribution, you "pool" the cells between both samples together, and then you randomly segregate the cells back into the two conditions maintaining original sample sizes. You then recalculate the proportional difference between the two conditions for each cluster, and compare that to the observed proportional difference for each cluster. I tend to take the log2 difference in proportions since it's a more sensible scale. Repeat this process about 10,000 times, and the p-value would be the number of simulations where the simulated proportional difference was as or more extreme than observed (plus one) over the total number of simulations (plus one).

Since I found myself having to do this so many times, I made a little R library for myself that takes a seurat object, and will do a permutation test for p-values (and adjusted p-values), as well as generate a plot with the observed proportional difference and a bootstrapped confidence interval for each cluster.

https://github.com/rpolicastro/scProportionTest
Question: Hi,

I would like to analyse the content of a GTF file. I am quite able with R and dplyr, so I would like to transform my GTF file into a data frame to facilitate my analysis. Does anybody know of any tool to do this?

Thanks.
Best, C.

Answer: Hi,

This worked:

gtf2 <- read.table('celegans.gtf', header = FALSE, sep = '\t')

gff <- read.delim('wormbase.gff3', header = FALSE, sep = '\t', skip = 8)

I forgot to specify the tab delimiter in the read.table() function. I thought that it was the default but it isn't.

The answer of cpad0112 is much better though because with his way, all the meta information of the 9th column is put in separate columns whereas with my way, all the meta information is all in one column.

Best,
C.
Question: Hi,

I have paired end data which im trying to align using tophat:

tophat -r 200 -p 4 -o /.../Tophat_output/ -G /pathtogtf /pathtobowtieindexfiles Sample10_untrim1.fastq Sample10_untrim2.fastq

Im getting only one bam alignment file. Im not sure if I should be getting 2 bam files.

Answer: You should get the accepted bam file, and another one of unmapped sequences

But don't expect two bam accepted files, one for each of the paired sequences. You only get one
Question: Hi! I am conducting a GWAS. I used plink 1.9 (--pca) to obtain first 10 population principal components which I entered in the regression model as covariates to correct for population stratification. I would like to the use the .eigenvec file, which is the output from plink --pca containing population principal components, to summarize the population structure and look for potential heterogeneity. I would be super grateful if anyone could help or refer me to a practical tutorial. Thanks a lot!

Answer: You know for sure that you need to include all 10 principal components as covariates? One should only include the principal components that are actually segregating your groups of interests and that are therefore likely to affect the statistical inferences that you make from your data. This is 'adjusting' for population stratification.

With a `.eigenvec` file, one can easily generate a principal components bi-plot for any pairwise combination of PCs:

R

setwd("/YourDir/")
options(scipen=100, digits=3)

#Read in the eigenvectors
eigen <- data.frame(read.table("plink.eigenvec", header=FALSE, skip=0, sep=" "))
rownames(eigen) <- eigen[,2]
eigen <- eigen[,3:ncol(eigen)]

summary(eigen)

#Determine the proportion of variance of each component
proportionvariances <- ((apply(eigen, 1, sd)^2) / (sum(apply(eigen, 1, sd)^2)))*100

plot(eigen[,1], eigen[,2])

legend("topleft", bty="n",
cex=1.5, title="",
c("Population 1","Population 2","Population 3","Population 4","Population 5"),
fill = c("yellow","forestgreen","grey","royalblue","black"))

#Sample image
[*need to manage colours and layout yourself*]

*[from: https://www.biostars.org/p/271694/]*

Kevin
Question: I've generated a volcano plot using DeSeq2 results and would like to specifically highlight a subset of genes by providing a list of gene IDs

Dataset$condition <- relevel(Dataset$condition, "Ctrl")
res <- lfcShrink(DatasetProcessed, contrast=c("condition","Treat","Ctrl"))
with(res, plot(log2FoldChange, -log10(pvalue), pch=16, cex=0.6, main="Volcano plot", xlim=c(-5,5)))

Given a list of geneIDs (corresponding to IDs in res) that I have converted into a vector, how can I then highlight these in the volcano plot?

gene_list <- scan("genelist.txt", what="", sep="\n")
gene_list <- as.vector(gene_list)

Answer: As always seems to be the case, I figured out the answer five minutes after asking the question...

with(res, plot(log2FoldChange, -log10(pvalue), pch=16, cex=0.6, main="Volcano plot", xlim=c(-5,5)))
gene_list <- scan("genelist.txt", what="", sep="\n")
gene_list <- as.vector(gene_list)
with(subset(res, rownames(res) %in% gene_list), points(log2FoldChange, -log10(pvalue), pch=20, col="orange"))

Hopefully this will help someone else down the line
Question: I've been using this method to rename files at the end of a step:

processed_fastq_1:
type: File
outputBinding:
glob: ${ return '**/' + inputs.fastq1.basename }
outputEval: |
${
self[0].basename = inputs.add_rg_SM + '_R1.fastq.gz';
return self[0]
}

But I'm not sure if this is an acceptable method. Should there be an intermediate step using python or calling to a script to rename files in the middle of a CWL workflow?

Answer: Hello ionox0,

That is acceptable way to rename, yes. However, `**` is not part of the POSIX glob definition, so you'll need to amend that. Perhaps the following simplification works with your tool?

processed_fastq_1:
type: File
outputBinding:
glob: */$(inputs.fastq1.basename)
outputEval: |
${
self[0].basename = inputs.add_rg_SM + '_R1.fastq.gz';
return self[0]
}

It has been requested that there be an easier way to rename output `File`s but there is no proposed new syntax yet. Perhaps you have a suggestion?

https://github.com/common-workflow-language/common-workflow-language/issues/668

The same link shows how to rename a `File` between steps.
Question: Hi,

I have some data that looks like this in R:

## ID td yr mo dy
## 1 1 0 18 3 2
## 2 1 1 17 3 12
## 3 1 2 16 1 4
## 4 2 2 18 3 7
## 5 2 0 19 2 8

I have some duplicates in column ID (there are sometimes triplicates) and would like the output to have those collapsed and to create new columns with the data being collapsed:

Would like:

## ID td yr mo dy td yr mo dy td yr mo dy
## 1 1 0 18 3 2 1 17 3 12 2 16 1 4
## 2 2 2 18 3 7 0 19 2 8

I have seen some solutions where they concatenate the data into a new column but would like to keep the data the same as above.

Thanks

Answer: Here is the start, using base **R**:

# example data
df1 <- read.table(text = "ID td yr mo dy
1 1 0 18 3 2
2 1 1 17 3 12
3 1 2 16 1 4
4 2 2 18 3 7
5 2 0 19 2 8", header = TRUE)

# split by ID, convert to vector
res <- lapply(split(df1, df1$ID), function(i) c(t(i[ -1 ])))

# set same lengths, then rbind
res <- do.call(rbind, lapply(res, `length<-`, max(lengths(res))))

res
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
# 1 0 18 3 2 1 17 3 12 2 16 1 4
# 2 2 18 3 7 0 19 2 8 NA NA NA NA
Question: Hi there,

I'm trying to create a pheatmap with annotationcolors based on groups. I'm an sociologist, so this might not be the common way to use pheatmap, but I found it to be great at displaying social groups as well. But,, when I try to plot an assymmetrical matrix, where the rows and columns do not match, I get an error:

```r
Error in check.length("fill") :
'gpar' element 'fill' must not be length 0
```

The issue is also tackled [here][1], but as I understand it, the solution is to make sure rownames and colnames are identical. Which isn't possible for my purpose.

here is an reproducible example:

```r
DAT <- matrix(c(
771,1176,1018,1042,275,949,823,1250,348,622,1244,
465,984,766,790,573,784,592,1250,1222,421,1103,
171,288,194,80,547,528,249,959,72,81,224,
132,122,90,50,44,206,127,108,52,59,144,
0,9,5,9,0,10,17,0,0,0,0),
nrow = 5, ncol = 11, byrow = TRUE)

group <- c("grp1","grp2","grp3","grp4","grp5","grp6","grp7","grp8","grp9","grp9","grp9")

annotation_c <- as.data.frame(group)

colnames(DAT) <- paste(group," ", "no.",seq_len(11),sep="")
rownames(DAT) <- paste(rep(c("grp"),5),91:95,sep="_")

# annotation_r <- as.data.frame(rownames(DAT))

pheatmap::pheatmap(DAT, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(seq_len(10))), cluster_cols=FALSE,cluster_rows=FALSE ,
annotation_col=annotation_c)
```

Thank you for your time,
Emil

[1]: https://www.biostars.org/p/175858/

Answer: I added rownames to the annotation_c data frame

rownames(annotation_c) <- colnames(DAT)

and your script worked.
Question: Hi,

In my abyss-pe command line, I use "np=10" since I am running abyss-mpi. At the mapping stag in which abyss-pe calls abyss-map, it work only on on thread!

Here is the leading part of the abyss-pe command:

$ /usr/local/abyss/1.9.0/bin/abyss-pe np=10 k=64 name=AbC -d \

Here is the abyss-map command that abyss-pe executes:

```
$ abyss-map -j1 -l64 P_R2_001_val_2.fq.gz P_R1_001_val_1.fq.gz Ab10_10-3.fa \
| abyss-fixmate -l64 -h pe150_1-3.hist \
| sort -snk3 -k4 \
| DistanceEst -j1 -k64 -l64 -s200 -n10 -o pe150_1-3.dist pe150_1-3.hist
```

Is there any way to make abyss-pe run abyss-map multi-threaded, i.e j10 for instance? I thought that abyss-pe will always make `j=np`, but this is not true. I am using abyss v1.9.

Thanks

Answer: Hi malabady,

Use `j=10` in addition to `np=10` in your abyss-pe command line.

`j` specifies the number of threads to use, for programs in the ABySS pipeline that are multithreaded (e.g. `abyss-map`). `np` specifies the number of MPI processes (a.k.a. "MPI ranks") to `ABYSS-P`, which is the first stage of the ABySS pipeline. `ABYSS-P` is the only program in the assembly pipeline that uses MPI and it is not multithreaded. The remainder of the pipeline runs on only a single machine.

Edit: fix typo and clarify that `ABYSS-P` is not multithreaded.

Edit 2: Correction: Use `j=10` *in addition to* `np=10`. Mention that the later stages of the pipeline only run on a single cluster node.
Question: Find a gene of interest in a species. I have a gene of a model animal for example (a certain gene in zebrafish), and I want to get this coding region (cds) of all fish for molecular evolution analysis. Due to the presence of introns, blasting the genome results in multiple hits (blast hits per exon). I need to manually adjust each hits and get the cds area I want, I want to ask if there is an easy way to achieve what I want.
Through a cds region, the corresponding cds region is obtained in the genome of the species of interest.![enter image description here][1]

[1]: /media/images/1a429411-4ee9-42c4-8632-fb1c9ab2

Answer: Easy is mostly relative to your skills, it is not trivial to do right though. Rule of thumb: the more conserved the protein sequence, the easier it becomes. For high sequence identities e.g. over 50% any approach is likely to work fine.

Easy: find orthologous sequences in fish species by using web-blast and extract the sequences of blast hits (given defined cut-offs) directly from GenBank. These will often be from predicted protein-coding sequences, which is fine if the annotation is at least half-reliable. This is all possible via the web-interface.

Easy: Extract orthologues via Ensembl's compara interface (gene tree) for the gene of interest. That will give you some orthologues, but only those in Ensembl

Intermediate: get the CDS instead of protein sequences or do the same on the commandline

Harder: if you want to include badly or unannoated genomes, you can use exonerate in est2genome or protein2genome mode. That may work well for well-conserved proteins. You can use your blast results to restrict contigs to search because exonerate alignment will be much slower than a blast run.

The following command will get you the CDS sequence:

exonerate --model protein2genome -q query.faa -t target.fna --refine full --ryo ">%qi vs. %ti %td %m aln.length: %qal score: %s %%-ident: %pi %%-sim: %ps strand: %tS CDS sequence \n%tcs"

Hardest: Annotate the unannotated genomes from scratch using MAKER2, BRAKER, etc. but I do not recommend it. The result for any single protein may not be better than the previous alternative.
Question: Can someone please help me count the number of reads in fastq file with 3' terminal bases? For example, I need to get the account of reads ending with A, T, G and C in any input.fastq

Answer: I like the `seqtk` version. You can also use a combination of native commands like `sed`, `rev`, and `cut` to accomplish the same thing and it's easier on the eyes.

sed -n '2~4p' input.fq | rev | cut -c1 | sort | uniq -c

`sort | uniq -c` is decent for low millions of lines. If your file is large, `awk` or `mawk` is generally preferred. You can replace it with `awk '{ cts[$0] += 1 } END { for (v in cts) print cts[v], v }'` and it'll save you quite a bit of time. If your file is really large, the `seqtk` version should be the fastest.

To use `awk` or `mawk` for sort unique, you can save it as an alias.

alias sortuniq="awk '{ cts[\$0] += 1 } END { for (v in cts) print cts[v], v }'"
Question: Can anyone tell me how to generate `vcf.gz` file and its index file `vcf.gz.tbi` in 1000 Genomes Project?

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20110521/

Answer: I have a script which does this using a VCF stream on stdin:

#!/bin/bash

file=$1

bgzip >$file
tabix -f -p vcf $file

I found I was always writing the same lines over and over when indexing VCF files. You can use it like this:

cat uncompressed.vcf | bgziptabix compressed.vcf.gz
Question: Assume we identify - by RNA-seq, tiling arrays, by prediction - possible candidate regions for non-coding, small RNAs. I wish to verify and predict the function of as many RNAs as possible by computation before going to the lab. One could use eg. Rfam to find similar sequences, after that we are left with more than 90% that have no match. One could predict the 2D structure using eg. [RNAfold][1], compare that using [RNAforester][2]. But that does not get me even close to a function prediction. Do you have experience with other tools or a better computational pipeline that gets more information out of the ncRNA candidates, possibly even something specific to bacteria.

[1]: http://www.tbi.univie.ac.at/~ivo/RNA/RNAfold.html
[2]: http://bibiserv.techfak.uni-bielefeld.de/rnaforester/

Answer: [?]

http://www.bioinformatics.fr/resources.php?tag=non-coding-rnas