text:
How to extract values from second file on the basis common first column?
File1:
ENST00000000233 ARF5
ENST00000000412 M6PR
ENST00000001008 FKBP4
ENST00000001146 CYP26B1
ENST00000002501 DBNDD1
ENST00000002596 HS3ST1

File 2:
ENST00000000412 0.779142
ENST00000001008 0.738143
ENST00000002501 0.715315
ENST00000002596 0.713664

Desired Output
ENST00000000412 0.779142    M6PR
ENST00000001008 0.738143    FKBP4
ENST00000002501 0.715315    DBNDD1
ENST00000002596 0.713664    HS3ST1

I like to use the join command for this. For your example above, you can simply use the following:
$ join file2 file1
ENST00000000412 0.779142 M6PR
ENST00000001008 0.738143 FKBP4
ENST00000002501 0.715315 DBNDD1
ENST00000002596 0.713664 HS3ST1

This assumes your files are sorted. If they aren't sorted, you can use input redirection if your shell supports it (bash and other modern shells do):
join <(sort -k 1b,1 file2) <(sort -k 1b,1 file1)

If your shell doesn't support <(), sort the files first:
sort file1 > file1.sorted; sort file2 > file2.sorted
join <(sort -k 1b,1 file2.sorted) <(sort -k 1b,1 file1.sorted)

But the number of lines in file 2 and file 1 are not same. In an output file, I need a 1st and 2nd column of file 2 and 3rd column values to be extracted from file1

The result matches your desired output - can you clarify what you want?

Note that the join command does not simply paste files together, it matches two files based on a common column (by default column 1), and then tacks on the additional columns from the files in the specified order. I think it is exactly what you want, given what you've written. One caveat is that the files must be sorted by the matching column, though in your example above, they already are.

If they're not sorted, this might work: join <(sort -k 1b,1 file1.txt) <(sort -k 1b,1 file2.txt)

@Scot Thank you for the answer, I figured out my file (completed one) was not sorted that's why it didn't work initially

Using awk, this can be done in the following way:
awk 'NR==FNR{a[$1]=$2; next}$1 in a{print $1, $2, a[$1]}' file1.txt file2.txt > same.txt

Explanation:
NR==FNR (along with the next command)

Will go through the first file only and create an associative array of of the geneIDs corresponding with the geneNames.
Then if the 1st column in the 2nd file is present in the array, a, it will print both columns from the file plus the geneName.


text:
Mixcr batch processing
I have a folder with 150 pairs of fastq files from an illumina sequence run. I need to precess them with mixcr, how can I do this in Bash with a single command?

What is your command for processing a single pair of FASTQ files?

You can use a for loop.

Here's an example of the command line that I'm using to run Mixcr for a pair of samples:  mixcr analyze amplicon --s hs --staring-material dna --5-end v-primers --3-end j-primers --adapters no-adapters --receptor-type trb --only-productive /Documents/Data/Input/S1_S1_L001_R1_001.fastq.gz /Documents/Data/Input/S1_S1_L001_R2_001.fastq.gz
The S1_S1_ component of the file name changes with each sample to S2_S2_  S3_S3_  all the way to 75.  I'm having a little problem setting the $variable for the file name in the loop.

One solution I thought is to try to iterate through the two types of files (R1 and R2) with a FOR loop, but I can't figure out how to do this.  I don't think a nested loop will do the trick.

Create a text file with two columns, one for R1 and the other for R2. Loop through this file. Shell loops are not a good thing to begin with, and nested loops are almost never required.

Here's how I'd approach this:
for READ1 in /Documents/Data/Input/*R1*
do
    READ2=${READ1/R1/R2}
    mixcr analyze amplicon --s hs \
        --starting-material dna \
        --5-end v-primers \
        --3-end j-primers \
        --adapters no-adapters \
        --receptor-type trb \
        --only-productive "$READ1" "$READ2"
done

Basically, loop over all the R1 FASTQ files in a directory, and calculate the R2 names as needed.
This is using BASH string manipulation, and assumes that substituting R2 for R1 will correctly calculate the R2 filename from R1. It also assumes that all the FASTQ files are in a single directory.

Thanks.  This is essentially the approach that I took.  My original problem actually was the result of trying to run a shell command that wasn't compatible with the Linux version of shell.  When I used $ batch instead of $ sh to initiate my script, all the recommendations worked. I should have started by providing the bit of information that I am working in Linux (Ubuntu 18.04 to be exact).

Ah - OK. Yes, you may want to use bash rather than sh if possible. Glad it's working for you.

I would also suggest using GNU parallel as described in our guides :
https://docs.milaboratories.com/mixcr/guides/generic-multiplex-bcr/#one-command-solution
Briefly:
ls /Documents/Data/Input/*_R1* | \
    parallel -j 4 \
    'mixcr analyze
    …
    {} \
    {=s:R1:R2:=} \
{=s:.*/:/Documents/Data/result/:;s:_R.*::=}'

That will take all fastq pairs from /Documents/Data/Input and put the result output files into /Documents/Data/result folder.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
extracting last few lines of a growing bgzip compressed .vcf.gz file
I am running bcftools merge to create a single .vcf.gz file out of a collection of chromosomal .vcf.gz files from the 1000g project.
The merge process has been running for a couple of days, and I want to estimate how much longer it will take. The file is a growing .vcf.gz that is being bgzip -c compressed from a Linux pipe:
bcftools merge [...] | some processing and filtering | bgzip -c > final.file.vcf.gz

How can I know at which point is the file currently, to predict how much longer it will take to complete the merge?
Is there a way to do something like:
tail --with-some-magic final.file.vcf.gz | bgzip -d -c

So that I can read at which chromosome is currently on?
Is there a way to guesstimate the final size of the merge file and compare from the sum of sizes of the original files? E.g. the sum of compressed sizes originally (before filtering and munging) is 20G:
find $PWD -name "ALL.chr*.vcf.gz" | xargs ls -l | awk '{print $5}' | csvtk summary -H -f 1:sum | numfmt --to-si

This gives 20G, but the merge file is at 27G and still running.
Any ideas?

27GB is not huge. You can zcat the file and see the last few lines. It will take ~10 minutes. It is theoretically possible to peek a bgzf file, but no such tools are available.
By the way,

I am running bcftools merge to create a single .vcf.gz file out of a collection of chromosomal .vcf.gz files from the 1000g project.

Probably you should use bcftools concat, not bcftools merge. merge combines VCFs from different samples. concat concatenates VCFs for individual chromosomes; each input VCF contains the same set of samples.

Ah! I didn't know about bcftools concat, I will definitely try that. The bcftools merge has been running for several days and no idea what's taking that long.

You could find which Chromosome the merged file is currently writing on by:
bcftools query -f'%CHROM\n' final.file.vcf.gz | tail

Note that also you might want to try multithreading the command that you are using to write the file, as bcftools is capable of using several threads for the compression of output files. E.g.:
bcftools merge -Ob --threads 32 file1.vcf.gz file2.vcf.gz > mergedFile.vcf.gz

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Where to find secondary metabolites in lichen (and herbal medicine)?
I've searched around and found about 80 specific secondary metabolites for lichen, but I can't seem to find anymore. They were mostly from research papers. I read that there are more than 800 secondary metabolites of lichen. Is there some website that has all of them? I tried using HMDB but it yielded nothing.

It is a hard question to answer. Could you provide some search terms you used to find the database. Also, please clarify if you want a database about metabolites of lichen or from herbal medicine. Also if you could please post the reference of the 800 secondary metabolites it would be helpful too to see how they got that number

I suppose that will also depend on the lichen species (complex). This would be perhaps more suited for biology.SE.

I'm voting to close this question as off-topic because it doesn't seem to be about bioinformatics at all.


text:
What do the FASTQ file names mean here?
I would like to run an analysis on Nanopore data. I'm trying to download sample data from https://github.com/nanopore-wgs-consortium/chm13. I would like to run minimap2 on the FASTQ files, resulting in a single BAM file.

Q: What do the file names mean? What is "S13"? What is "L002_l1"? What is "L002_R1"? Am I supposed to download them and combine them into a single FASTQ file, and then align them?

What is "S13"?
S# is most likely the sample number. If so, you do not want to combine these.
What is "L002_l1"?
L00# is the lane number (of the 8 Illumina flow cell lanes). 'I1_' (not 'L') are index files and are not needed (see also https://www.biostars.org/p/374581/).
What is "L002_R1"?
R1 is one end of paired-end sequencing, R2 is the other end. These filepairs must have the same number of lines. Use an aligner that utilizes both.

Peter's answer is a good one and to the point. I'll be a bit more verbose.
The breakdown of FASTQ file names that come directly from the sequencer typically have the following format:
{sample_name}_S{sample number}_L{lane number}_{R/I}{read or index number}_001.fastq.gz
sample_name is the sample name provided by you (or whoever sequenced the data) to the sequencer.
To save time and money, samples are often "multiplexed" on sequencing runs. This multiplexing means you aggregate samples ("S") on the same run and split samples on different lanes ("L") to limit some technical noise. Lanes are essentially technical replicates for the sequencing data on the same sequencing run, and these numbers tell you which sample and which lane that FASTQ comes from.
Sequencing can be done in single-end or paired-end. Single-end sequencing generates a single FASTQ file (_R1), whereas paired-end sequencing produces 2 FASTQ files (_R1 and _R2).
The data you highlight isn't actually Nanopore:

Approximately 50x of data was generated on a NovaSeq instrument.

10x includes "indexes" ("I") in addition to the sample and lane information to keep track of things like single cells or other unique molecular identifiers, and that information is listed after the I. _R1, _R2, and _I1 all go together for a unique sample and lane combination

The file names mean this isn't Nanopore data. Read the bit just above these file names:

Approximately 50x of data was generated on a NovaSeq instrument.
Based on the summary output of Supernova, there are 1.2 billion reads with 41x effective coverage.

That linked github page has Nanopore, Illumina/10x, BioNano, HiC, and PacBio data. For the nanopore files, you should be looking at this section.


text:
How to extract all pair end fastq.gz files from multiple subdirectory to a single folder
I have all file names for each paired-end fastq pair as shown below. Could you please suggest how I can extract all fastq.gz files in a single folder?
sample name:
188833984]$ ls

Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc
Leaf_T1_FD_R2_L001-ds.7db8eb7e3426486db549426601b3a0bd
Leaf_T1_FD_R3_L001-ds.147177ecc03a46ccbbce3162d1185a0a

pair end file within sample name Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc
188833984]$ cd Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc

Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc]$ ls
Leaf-T1-FD-R10_S73_L001_R1_001.fastq.gz
Leaf-T1-FD-R10_S73_L001_R2_001.fastq.gz

I need all in fastq.gz format but the problem is that each sample, pair end R1 and R2 fastq reads in a single folder like this: and I need to extract all samples files in a single folder

Can you please [edit] your question and clarify. What do you mean by "extract all fastq files"? Do you just want to decompress (go from fastq.gz to fastq)? What is the final outcome you are expecting from the example you give?

I think OP wants to decompress fastq.gz files found in various subdirectories to fastq files in a single location. OP should probably decompress them where they are and create soft-links to them in their directory of choice.

@Dr.Yogesh, Do you want to copy/move all files to one folder or decompress them (as Ram RS is suggesting)?

"this file name has R1 and R1 fastq.gz": what does this mean? Did you mean to write R1 and R2?  You show a directory, Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc with two files, an R1 and an R2, so what is the problem? The files seem to already be in separate directories. Can you please [edit] and show us what the final result you want is?

Go to the common parent folder under which all FASTQ.GZ files reside. They could be in sub-folders under the folder, but all of them should have the folder you're in as a common parent folder. From such a folder, you should be able to run a find
find . -name "*.fastq.gz"

to get your FASTQ.GZ files. Now, there are surely better ways of doing this but here's one straightforward way.
Let's say you need to extract the FASTQs to ~/target/dir/, you'd do something like:
for f in $(find . -name "*.fastq.gz")
do
    tgt_file="~/target/dir/${f%.*}" #removed .gz from the name
    gzip -dc $f >$tgt_file #pigz -dc would be better if you have pigz
done

EDIT: OP has clarified that they wish to copy the files as such, not decompress them.
There are two ways to achieve this goal - you can either make physical copies of the files (thereby doubling the storage space you're currently using for them), or you can create symbolic links to the files (which are like shortcuts pointing to the actual files). I'll show you both ways here:
Let's assume the files are all under /fastq/files
find /fastq/files/ -name "*.fastq.gz" -exec cp {} ~/target/dir/ \; #This creates physical copies

cd ~/target/dir
find /fastq/files/ -name "*.fastq.gz" | xargs -I v_f ln -s v_f #These lines create symbolic links that point to the original files. Do not delete/move/rename the original files or these will become invalid.

Dear Ram, I want all subdirectory fastq.gz files into a single folder but files should be in fastq.gz format only, I don't want to decompress the files.

I've updated my answer. Use symbolic links instead of copying the files. As long as you don't move the original files or rename them (or any of their parent directories), your links will work fine and take much less space than if you copied the actual files.

@Dr.Yogesh Please accept the answer if it worked for you.


text:
Run Rscript for multiple samples (files)
I have used the following script for calculating the median for one file (gene). However, I have 20,000 files for which I want to run this R script and get 200 separate output files
 #!/usr/bin/R
 setwd("/location folder")

 d <- read.table("tgf.txt", header=FALSE)

 x <- aggregate(d[, 2], list(d$V3), median)

 write.table(x, "tgf.output.txt", quote=FALSE, row.names=TRUE)

 q()

Another option is to use list.files and grep to select only the .txt in your directory. A for loop can read-aggregate and save the files one at the time.
all.files <- list.files("/location folder")
txt.files <- grep(".txt",all.files,value=T)

for(i in txt.files){
  d <- read.table(paste("/location folder",i,sep=""),header=F)
  x <- aggregate(d[, 2], list(d$V3), median)
  write.table(x, paste0(i,".output.txt"), quote=FALSE, row.names=TRUE)
}

Why not just list.files(pattern = "*.txt")? The grep is not required.

And looping when vectorized operations are available is de-optimization.

And there's some inconsistent usage: paste(..., sep="") on line 5 and paste0() on line 7.

considering the nature of the question I thought that it was more important to show possibilities rather than one single point of view. paste and paste0 achieve the same but in different ways (e.g. sep="").  I tend to split commands involving list.files and grep as it makes easier to read and to reuse. The content of all.files is hard to miss. Finally the usage of the applyfamily is not easy to grasp for new user, while loops are much intuitive. Furthermore, apply is not always better than a nice for, especially in an example like this.

Showing possibilities is great! However, I don't think a loop would perform as well as or better than vectorized operations. While a loop is easier for people that know basic programming, the apply family solves a lot of problems, especially when there is a learning opportunity such as this one, where the task is well defined and how apply works can be easily understood.

Just encountered this particular gotcha, so we might have to be careful about using T/F instead of TRUE/FALSE as well: https://www.johndcook.com/blog/r_language_for_programmers/#assignment

Please do not use setwd in a script, it changes working folder without explicit specifications from the user, so the results could be unpredictable for people that are not familiar with the content of the script.
Create a tab-separated file with input filename and output filename, or at least a file with all input files, unless you can pick them using a list.files(pattern = "<pattern_here>"). Read this file as a data.frame and use it with an apply, so a function processes the df row-by-row. Let's say the file (named metadata.txt) is:
in_file out_file
tgf.txt tgf.output.txt
abc.txt abc.filtered.txt
xyz.txt xyz.randomname.txt

Define a function called process_row() in the Rscript:
process_row <- function(r) {
    d <- read.table(r[[1]], header = FALSE);
    x <- aggregate(d[, 2], list(d$V3), median);
    write.table(x, r[[2]], quote=FALSE, row.names=TRUE);
}

And do the actual work:
meta_data <- read.table('./metadata.txt', sep = ' ', header = TRUE, stringsAsFactors = FALSE);
apply(meta_data, 1, process_row(r))
# You don't need a q() in an Rscript because it exits after the last line of code anyway.

I tried the above script but it is giving me following error
Error in read.table(r[[1]], header = FALSE) :
'file' must be a character string or connection
Called from: read.table(r[[1]], header = FALSE)

What is the content of ./metadata.txt, which is in the working directory where you're running R?

It worked. Thank you

Please accept the answer if it worked. Thank you


text:
How can I translate infix to mathml using the python interface to libsbml?
Assume I have an infix expression such as k1*s1, how do I use libsbml to translate this expression to mathml using the python interface?

Welcome to the site.  Can you include examples if what you have already tried?

See answer below. It works, I tried it.

It’s quite straightforward to convert infix to mathml using libsbml. Make sure you have libsbml installed in python by typing
pip install python-libsbml

at the windows/Mac/Linux command line or
!pip install python-libsbml

if you are running a Jupyter notebook. Some tools such as Tellurium already have libsbml preinstalled.
At the python command line type
import libsbml

to import the package. To enter the infix expression type the following:
p = libsbml.parseFormula ('k1*s1')

To generate the mathml type
print (libsbml.writeMathMLToString (p))

One should use the newer libsbml.parseL3Formula which provides the full set of SBML Level 3 formula parsing.  I recommend using the extended libsbml.parseL3FormulaWithModel(formula, model) which uses a model instance to resolve the actual symbols in the formula. Thereby it is clear if all the symbols exist in the model. If the functions return None make sure to check the getLastParseL3Error to find out what went wrong.


text:
UCSC liftover cannot produce bed format output
I was trying to use UCSC liftover tool to convert human coordinates to mouse. While most of the files worked, one file didn't:
chr16   86430087    86430726    Human_chr16:86430087-86430726___element_1___positive____neural_tube[12/12]___hindbrain_(rhombencephalon)[12/12]___limb[3/12]___cranial_nerve[8/12]
chr16   85620095    85621736    Human_chr16:85620095-85621736___element_2___negative_
chr16   80423343    80424652    Human_chr16:80423343-80424652___element_3___negative_
chr16   80372593    80373755    Human_chr16:80372593-80373755___element_4___positive____neural_tube[6/10]___hindbrain_(rhombencephalon)[10/10]___midbrain_(mesencephalon)[10/10]

Does anyone know why this bed file can only be converted to chr:start-end format but not bed4 format?

Were there any logged messages during the process?  Can you include an example that did complete the liftover?

If you remove everything after the "Human" in each line then it works

I figured out the issue by myself. It turned out that UCSC liftover doesn't like the character ":"
I removed ":" and it worked


text:
Subsetting of object existing of two samples
I have a Seurat object existing of an aggregate of two samples namely; RD1 and RD2. I am trying to make a subset of each sample. But it's not really working. I manage to make a subset of the cells but I don't manage to make a subset of the whole object. Isn't there a way to do this, based on the subset from cell names a subset to make from the whole object.
seuratObj$orignal_1 <- seuratObj@meta.data$orig.ident
seuratObj$orignal_2 <- seuratObj@meta.data$orig.ident
seuratObj$orignal_1 <- gsub("RD1", "TRUE",seuratObj$orignal_1 )
seuratObj$orignal_1  <- gsub("RD2", "False",seuratObj$orignal_1  )
seuratObj$orignal_2  <- gsub("RD1k", "TRUE",seuratObj$orignal_2  )
seuratObj$orignal_2  <- gsub("RD2", "False",seuratObj$orignal_2  )

subset_DR1 <- data.frame(orignal_1 = seuratObj$orignal_1)
subset_DR1 <- subset(subset_DR1, orignal_1 == TRUE,select = c("orignal_1"))

subset_DR2 <- data.frame(orignal_2 = seuratObj$orignal_2)
subset_DR2 <- subset(subset_DR2 , orignal_2 == TRUE,select = c("orignal_2"))

SubsetData(seuratObj, cells = rownames(subset_DR1), subset.name = seuratObj$orignal_1)

I just started to work with seurat packages. I think I am missing something.

Welcome to the site ageans. Did you try to create two seurat objects one for each sample? (Instead of dividing this single seurat object?)

Yes, I have 1 suerat object and now I am trying to divide them into two Seurat object, 1 for each sample.

Ok, but can't you create from the start two seurat objects?

Sadly enough now because it's an old object. And I only got the object. And for the rest of the analysis, it should be 1 object. But for 1 part I need to separate it.

Can you post some information on the meta data? Do you have a column in the metadata that annotates the subset?Also assign the subset to new object, subset_obj <- SubsetData() if you use seurat2, subset_obj <-  subset(x = seuratObj, subset = orig.ident == "Replicate1") if you use seurat3.

I think your issue is simply not assigning the newly-subset object in to a new object.
Assuming rownames(subset_DR1) is indeed a list of cell names
seuratObj_subset_dr1 <- SubsetData(seuratObj, cells=rownames(subset_DR1))

will work.


text:
How to do `bedtools intersection` using pandas alone?
I have two pandas Dataframes, using python3.x:
import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)
print(df2)
##       0    1    2         3  4  5
## 0  chr1  155  200  feature5  0  -
## 1  chr1  800  901  feature6  0  +

The columns to focus on in these dataframes are the first three columns: location, start, and end. Each start:end value represents a distance on location (e.g. chr1, chr2, chr3).
I would like to output the intersection of df1 against df2. Here is the correct output:
chr1    155 200 feature2    0   +
chr1    155 200 feature3    0   -
chr1    900 901 feature4    0   +

What is the most efficient (in terms of runtime and RAM) to find the intersections using pandas?
There are python wrappers around BEDTools: https://daler.github.io/pybedtools/intersections.html

I was hoping there was a pandas only solution...

https://github.com/biocore-ntnu/pyranges is likely what you want. I'll write an answer :)

Your setup:
import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)

Now use pyranges (https://github.com/biocore-ntnu/pyranges):
# pip install pyranges
# conda install -c bioconda pyranges

import pyranges as pr
cols = "Chromosome Start End Name Score Strand".split()
# Add column names to the dfs so PyRanges knows which
# represent Chromosomes, Starts, Ends and optionally Strand

df1.columns = cols
df2.columns = cols

# create PyRanges-objects from the dfs
gr1, gr2 = pr.PyRanges(df1), pr.PyRanges(df2)

# intersect the two
gr = gr1.intersect(gr2)

# get the result as a   df
df = gr.df

Results:
print(gr)
# +--------------+-----------+-----------+------------+-----------+--------------+
# | Chromosome   |     Start |       End | Name       |     Score | Strand       |
# | (category)   |   (int32) |   (int32) | (object)   |   (int64) | (category)   |
# |--------------+-----------+-----------+------------+-----------+--------------|
# | chr1         |       155 |       200 | feature2   |         0 | +            |
# | chr1         |       900 |       901 | feature4   |         0 | +            |
# | chr1         |       155 |       200 | feature3   |         0 | -            |
# +--------------+-----------+-----------+------------+-----------+--------------+
# Stranded PyRanges object has 3 rows and 6 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.
print(df)
#   Chromosome  Start  End      Name  Score Strand
# 0       chr1    155  200  feature2      0      +
# 1       chr1    900  901  feature4      0      +
# 2       chr1    155  200  feature3      0      -

See the docs for more info: https://biocore-ntnu.github.io/pyranges/

Hi @TheUnfunCat , looks very promising. I haven't try out the full package. I look forward to reading the publication. In your example, ever thought of including what are the intersecting features from gr2? In this case should be feature5, feature6, feature5 ..

Thanks for this answer! I feel this is a really useful contribution, so I marked it as correct---naturally, others can ask as well. pyranges looks great!
One request: to improve your answer, could you add comments? That would be helpful for new readers unfamiliar with the pyranges package. Thanks!

@StupidWolf Use gr1.join(gr2) for that :)

pyranges looks good, I see your benchmarks have it orders of magnitude faster than GenomicRanges in R.  I almost want to try python instead of using GenomicRanges in R, almost.

Agreed pyranges is a brillant solution by @TheUnfunCat

I don't think Pandas has this implemented functionality out-of-the-box. Even if it did, solutions not designed specifically for bioinformatics probably rarely handle intervals on different chromosomes correctly unless you split the intervals by chromosome first.
Pandas does handle intervals (see docs for the Interval and IntervalIndex classes), but I've never used these. I have used the excellent intervaltree package. Whatever approach and tools you decide to use, there's a good chance you'll have to write some code. Computing interval intersections with intervaltree would look something like this.
>>> from collections import defaultdict
>>> from intervaltree import IntervalTree
>>> from pandas import DataFrame
>>>
>>> dict1 = {
...     0:['chr1','chr1','chr1','chr1','chr2'],
...     1:[1, 100, 150, 900, 1],
...     2:[100, 200, 500, 950, 100],
...     3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
...     4:[0, 0, 0, 0, 0],
...     5:['+','+','-','+','+']
... }
>>> dict2 = {
...     0:['chr1','chr1'],
...     1:[155, 800],
...     2:[200, 901],
...     3:['feature5', 'feature6'],
...     4:[0, 0],
...     5:['-','+']
... }
>>> df1 = DataFrame(dict1)
>>> df2 = DataFrame(dict2)
>>>
>>>
>>> # Construct an interval index from one of the dataframes,
>>> #   one tree per chromosome
>>> index = defaultdict(IntervalTree)
>>> for n, row in df2.iterrows():
...     index[row[0]].addi(row[1], row[2], row)
...
>>>
>>> # Query the index from the other dataframe to build a 3rd
>>> #   dataframe representing the intersection
>>> dict3 = { i: list() for i in range(6) }
>>> for n, row in df1.iterrows():
...     overlapping = index[row[0]].overlap(row[1], row[2])
...     for itvl in overlapping:
...         itsct_start = max(row[1], itvl.begin)
...         itsct_end = min(row[2], itvl.end)
...         dict3[0].append(row[0])
...         dict3[1].append(itsct_start)
...         dict3[2].append(itsct_end)
...         dict3[3].append(row[3])
...         dict3[4].append(row[4])
...         dict3[5].append(row[5])
...
>>> df3 = DataFrame(dict3)
>>> df3
      0    1    2         3  4  5
0  chr1  155  200  feature2  0  +
1  chr1  155  200  feature3  0  -
2  chr1  900  901  feature4  0  +
>>>

The Intervaltree package is great for some things, but if you want something > 100 x faster that works directly with pandas data please look into my https://github.com/biocore-ntnu/ncls library :)

As mentioned by OP, another option is to use pybedtools, which in my opinion is pretty convenient for people already familiar with BedTools.
Let's even say df1's format is slightly different than df2:
import pandas as pd

dict1 = {0: ['chr1', 'chr1', 'chr1', 'chr1', 'chr2'], 1: [1, 100, 150, 900, 1],
         2: [100, 200, 500, 950, 100],
         3: ['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
         4: [0, 0, 0, 0, 0], 5: ['+', '+', '-', '+', '+'],
         6: ["remember", "the", "5th", "of", "november"]}

df1 = pd.DataFrame(dict1)

print(df1)

      0    1    2         3  4  5         6
0  chr1    1  100  feature1  0  +  remember
1  chr1  100  200  feature2  0  +       the
2  chr1  150  500  feature3  0  -       5th
3  chr1  900  950  feature4  0  +        of
4  chr2    1  100  feature4  0  +  November

dict2 = {0: ['chr1', 'chr1'], 1: [155, 800], 2: [200, 901],
         3: ['feature5', 'feature6'], 4: [0, 0], 5: ['-', '+']}

df2 = pd.DataFrame(dict2)

print(df2)

      0    1    2         3  4  5
0  chr1  155  200  feature5  0  -
1  chr1  800  901  feature6  0  +

Using pybedtools we get:
from pybedtools import BedTool

bed1 = BedTool.from_dataframe(df1)
bed2 = BedTool.from_dataframe(df2)

intersected_bed = bed1.intersect(bed2)

print(intersected_bed)

chr1    155 200 feature2    0   +   the
chr1    155 200 feature3    0   -   5th
chr1    900 901 feature4    0   +   of

We can then keep using the intersected BedTool object or convert it back into a DataFrame:
# create a pandas.DataFrame, passing args and kwargs to pandas.read_table
intersected_df = intersected_bed.to_dataframe(header=None,
                                              names=["chrom", "start", "end", "name",
                                                     "score", "strand", "whatever"])

print(intersected_df)

  chrom  start  end      name  score strand whatever
0  chr1    155  200  feature2      0      +      the
1  chr1    155  200  feature3      0      -      5th
2  chr1    900  901  feature4      0      +       of

Thanks for the answer! It's been awhile since I've tried pybedtools. Isn't it heavy with respect to RAM?

I'm definitely not an expert at this, so please correct me if I'm mistaken, but according to this performance test it seems that while their memory usage is pretty much the same, running time is where pyranges beats pybedtools.
@TheUnfunCat

I implemented pandas "Intervals" and  ... it should be a few lines, clearly there are limitations. For non-overlapping data it is very cool however. It will work for overlapping data, BUT if the data you are using as the interval data is overlapping, it falls over. It could work if an independent (non-overlapping) interval was constructed.
Anyway, the point of the excercise was to keep it all in pandas,  "start_stop" should not be needed. This can't be done via Intervals, you've got to drop out, and that will slow it down and harder to parse.
The query contig df3.iloc[1,], is still a bug (hence not printed), should be [x,] where x is the value_counts hit, but this is just a pilot.
The script would be fine if it the df was wrangled using groupby rather than value_counts because value_counts trashes important data, which is what I'm trying to recover via df3.iloc[1,]. I like pandas elegant solutions, but you know another day another dollar.
Note: I changed 200 to 199 otherwise there's too much overlap.
Overlap Output (not labeled very well)
                 start  stop
index
(1.0, 100.0]        0     0
(100.0, 199.0]      1     0
(150.0, 500.0]      1     1
(900.0, 950.0]      0     1

import pandas as pd

def tableit (dict):
    df = pd.DataFrame(dict)
    table = pd.pivot_table(df, values=['start', 'stop'], columns=['chromosome'], index=['feature', 'misc', 'misc2'], fill_value = '-')
    return table

def start_stop(point, table2, bin):
    df3 = pd.cut(table2[point,'chr1'], bin[0]).value_counts()
    df3 = df3.to_frame(point).reset_index()
    for x in bin:
            if not (x == bin[0]):
                df3_tmp = pd.cut(table2[point,'chr1'], x).value_counts()
                df3_tmp = df3_tmp.to_frame(point).reset_index()
                df3 = pd.concat([df3, df3_tmp]
            else:
                continue
    df3.set_index('index', inplace=True)
    return df3, df3.iloc[1,] # [1,] is a bug and the query locus needs correctly reading from df3_tmp

def main():
    dict1 = {'chromosome':['chr1','chr1','chr1','chr1','chr2'], 'start':[1, 100, 150, 900, 1], 'stop':[100, 199, 500, 950, 100], 'feature':['feature1', 'feature2', 'feature3', 'feature4', 'feature4'], 'misc':[0, 0, 0, 0, 0], 'misc2':['+','+','-','+','+']}
    dict2 = {'chromosome':['chr1','chr1'], 'start':[155, 800], 'stop':[200, 901], 'feature':['feature5', 'feature6'], 'misc':[0, 0], 'misc2':['-','+']}

    table1 = tableit(dict1)
    table2 = tableit(dict2)

    bin = [[i, j] for i, j in zip(table1['start', 'chr1'], table1['stop', 'chr1'])]

    df_start, locus = start_stop('start', table2, bin)
# locus.name gets to the data
    df_stop, _ = start_stop('stop', table2, bin)
    df_merge = pd.merge(df_start, df_stop, on="index", how='outer')
    print (df_merge)

if __name__ == "__main__":
    main()

The following would have been better,
table_join = pd.concat([table1, table2]).fillna('-')

and then sort them.. perhaps including a new column to specify the original dataframe and screen for the overlap with table_join.iloc[i, column] - table_join.iloc[1-i, other_column].
The key step in all approaches is to convert from "long format" to "wide format" and thats the most important thing I've done in all code and makes pandas work.
Output Table_join
                     start        stop
chromosome            chr1 chr2   chr1 chr2
feature  misc misc2
feature1 0    +        1.0    -  100.0    -
feature2 0    +      100.0    -  199.0    -
feature3 0    -      150.0    -  500.0    -
feature4 0    +      900.0    1  950.0  100
feature5 0    -      155.0    -  200.0    -
feature6 0    +      800.0    -  901.0    -

At some point in the future I'll finish the script (its easier than starting over in this instance). The issue is parsing "pd.cut(table2[point,'chr1'], x)", assigning a forth index, which is unique to all rows, and then reading the output against the new index of the original table2 dataframe to obtain the dessired result. The +1 nucleotide issue is a pain, but anyway. Its only approx. 4 lines of code.


text:
Looking for a taxonomic dataset that links tree genera (or species) to families
I am looking for a taxonomic dataset that would allow me to links tree genera (or species) to families. In particular, I am interested in the tropical species found in the Amazon.
I found a long but not exhaustive list on Wikipedia:
https://en.wikipedia.org/wiki/List_of_trees_and_shrubs_by_taxonomic_family
After processing it, I realised it contains information about 404 genera, while the dataset I am dealing with has around 600. (in the end, I only got 23% of my data covered, very sad)
Does anyone know if such a crosswalk file exists, in some scientific article or somewhere in the Internet?
Moreover, please tell me if there are other StackExchange websites which could be more appropriate, I wasn't sure of where to post this.

Contact Kew. BTW when emailing them make sure you call them genera.

Thanks, it is always good to learn something new!

Probably the best currently available online database is RBG Kew's Plants of the World Online. It organises taxonomic, morphological and distribution data from a range of databases. You can search for your taxa of interest and find their associated higher ranks.
Alternatively, there is World Flora Online, replacing the now-defunct The Plant List v1.1, which also integrates taxonomic data from a range of reliable sources and national lists.

Thanks! RBG Kew's site looks great, do you know if I can somehow download the data?

Hi @Davide - you should use the dedicated APIs: either pykew for python, or taxize for R.

Thanks @NatWH, I'll look pykew up

Hey @NatWH, I am having problems with using pykew, and I cannot find many usage examples. To do what I am asking I try  to use the function powo.search({Name.genus = 'Poa'}), for example, and I just get the answer <pykew.core.SearchResult object at 0x00000231EBA12DA0>, and I don`t know how to extract the family name out of it...

@Davide pykew returns a search result as an iterator consisting of as many dictionary-like elements as your search term has hits. Try this fragment to get the family associated with each record: res = powo.search({Name.genus: 'Poa'}) followed by [ r['family'] for r in res ]. Irritatingly (but understandably), powo.search() returns a generator rather than an iterator, so you will have to rerun the search if you want to operate on the same results again.

Thank you @NatWH, now I manage to get some results! While running, I received an error when using your code for the genus 'Zizyphus'. I get " 'SearchResult' object has no attribute '_results' ". I imagine that such genus is not in the Kew database - in this case how can I tell the code to ignore the error and just return something like 'Unassigned'?

Actually realised that my data is not perfect: the correct spelling of that genus is 'Ziziphus' - pykew works fine. Thanks again!

The NCBI taxonomy is a quite extensive resource: https://www.ncbi.nlm.nih.gov/taxonomy
I would thus expect there is a good chance that it covers all of the taxa that you are interested in.  It officially claims that it is not an authoritative source, but it's still pretty darned good.

Please note that you should take the NCBI seriously when they say it is not an authoritative resource. Orthography and synonymy are rarely up to date, and many accepted species are not present because they lack associated sequence data.


text:
Effect size in power analysis when dealing with Poisson based variant caller
I was trying to estimate the minimum sequencing depth required to confidently detect somatic variant using Poisson distribution based variant caller (such as somatic variant caller or Pisces from illumina).
Power analysis is used which requires significance level (0.05), power / sensitivity (99%) and effect size to estimate sample size (sequencing depth). I am not sure how to calculate the effect size in such context. Can anyone share some insight? Thanks a lot.

It's standard stats, it's a while since I've done this, but it's textbook stuff

Try running a pilot study to get some preliminary results you can use to estimate the effect size = difference / SD.

authors: ə̷̶̸͇̘̜́̍͗̂̄︣͟


text:
Differences in webtool and standalone Clustal Omega Percent Identity Matrix
I'm trying to align sequences using MUSCLE in my own computer. To compute Percent Identity Matrix I'm trying to use Clustal2.1 like it's done in https://www.ebi.ac.uk/Tools/msa/muscle/. Even though I'm doing the same thing, I get rounded up values in my computer.
For example, this is the webtool output:
#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1  100.00   92.79   65.29   63.98
     2: M97264.1     92.79  100.00   65.66   63.92
     3: AJ292230.1   65.29   65.66  100.00   66.88
     4: KM379144.1   63.98   63.92   66.88  100.00

And this is the output in my computer.
#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1     100      93      65      64
     2: M97264.1        93     100      66      64
     3: AJ292230.1      65      66     100      67
     4: KM379144.1      64      64      67     100

I have both installed clustalw using sudo apt install clustalw and by downloading the tarballs from http://www.clustal.org/download/current/clustalw-2.1.tar.gz, with exact results.
The commmand I'm using is clustalw2 alignment -tree -pim.
Any thoughts on why this could be? Thanks in advance.

So, after a year I came out with the solution - as M__ kindly suggested, using Clustal Omega.
First, align usign muscle:
muscle -in fasta -quiet -out alignment
Then, generate percent distance matrix using clustal omega:
clustalo --infile alignment --distmat-out output --percent-id --full

Try clustal omega, instead of clustalw because it is really old. It should report to 2 dp. Variations do occur depending on how indels are treated, but this is just rounding.
You could try
alignment -tree -pim -score

The idea is to convert to absolute rather than percentage and might overcome the rounding error


text:
Is there a way to use BioLemmatizer from Python?
Lemmatization is one of the most important tasks for text normalization in natural language processing (NLP), i.e. using deep learning to understand text data.
One of the good lemmatizers I found for biomedical domain is BioLemmatizer. However, it is developed in Java and is available from maven repo. I would like to use it from Python, but there seems to be no counterpart for BioLemmatizer in Python. Is there a way I can access BioLemmatizer from Python using some module like JPype. JPype, however, supports access only for the default classes in Java.

This area is huge in Python because of machine learning text processing, in fact O'Reilly have an entire book on it . However, we don't really deal with words in Bioinfo SE,
You might want to look at the nltk (natural language toolkit) for Python, but for direct parallels with the Java package bio package (your question) I don't think it's something we would know.

You could solve this via subprocessing into Java, and passing the result back to Python. I get the idea that the core ML/ DL used in natural language is set up for Python. Google for example have fantastic tensorflow stack for this.

Unfortunately, bio package has many tools dealing with biological sequences, but I require a lemmatizer that work with biomedical literature in English.

Hi I'm aware of your goals, lemmatization orietated towards biological literature. Its not something we're clued up about, but it's an advanced area of programming


text:
Normalization of single cell RNASeq data with ERCC spike-ins
I wish to normalize a scRNASeq dataset with respect to ERCC spike-ins, where "for some of the samples, ERCC spike-in RNA was added to the lysis buffer" and was wondering how to do so?
I have seen similar questions posted on the forum, such as the following: Normalization methods with RNA-Seq ERCC spike in? and have checked out the package RUVSeq in the top answer. However, in this package, it seems that spike-ins are taken as transcripts, whereas in the dataset I'm analyzing it seems as though 'spike-ins' are respective of the cell being sequenced rather than a transcript. In addition, in my research I stumbled upon the following tutorial for preprocessing of scRNASeq data: https://f1000research.com/articles/5-2122/v2 where also, as in the RUVSeq package, "ERCC" is grep'd through the rownames to identify ERCC transcripts.
Which leads me to my second question: what is a spike-in if it is not a transcript - and therefore not within the rownames of the dataset and rather within the samples? Please note that I am not from a biological background and, I have read the wikipedia article, along with many others, in my quest to understand what a spike-in is.
In addition, please note that I have a metadata file that states exactly which of the samples have these spike-ins and which dont.
Any advice on how to proceed would be much appreciated.
Thanks.

Please use the search function. Normalization methods with RNA-Seq ERCC spike in?

@story, If you read the full text you will be able to see that I have already referred to that post

You should treat your ERCC's as entries in your table of "gene counts". You should have counts for them as individual rows (one for each ERCC) in your expression table. The easiest way to do this is to include them in your alignment by appending them to the genomic fasta file (used in the alignment or index generation).

I do not have the raw files, just the final count matrix, with metadata identifying which samples/cells had ERCC spike ins added to the lysis buffer.. Is there an alternative way to do what you are referring to for such instances? I have found nothing..

You need to get a quantitative readout of the ERCC's. You could try counting them directly (without alignment) from the .fasta/.fastq files. But you need those original files. The whole idea of adding ERCC's is to control for variation in yield... identical amounts should have been added to each sample that contained them; using that information you can then gauge whether there were major differences in concentrations in your final samples. Are you sure that the count table you have doesn't include them? Seems like significant oversight for whoever provided those files.

OK I understand. And yes, I'm sure that the count table does not include them. I will message the authors and ask for that information. Thanks for your input btw!


text:
How to get a MSA fasta from BAM/SAM?
I am working with paired end NGS miseq data from a viral genome with multiple timepoints.  I have filtered and trimmed this data for quality and adapter sequences.  I have then merged the filtered paired end reads and mapped them to a reference genome using bowtie2.  I then used samtools to filter for mapping quality and visualized the resulting filtered bam file on IGV.
My goal is to now scan along the genome using various window sizes (e.g. 300 nt) across the reads that have mapped to the reference. I will then chop up the genome and build phylogenetic trees based on each window.
How can I extract a multiple sequence alignment of these 300 nt windows from a sam or bam file using the coordinates of a reference viral genome?  I have tried using samtools view but this does not actually provide the subsection or alignment within the window but rather all reads that overlap the window in some way.
The virus is HIV1 and the immediate goal is to get an alignment of all reads that completely span a user specified coordinate range eg. 1-300 based on some reference sequence. Ideally reads that fit this criteria would be extracted and aligned. After ragged edges are removed, the resulting alignments subjected to phylogenetic analyses.

this seems a little "chaotic". I would consider doing variant calling and then trying some approach?

If need to get msa between every read and ref in the specific region ,
you can try this https://github.com/orangeSi/bam2msa.git
bam2sam ref.fa input.bam --regions chr1:200-300


text:
Filter reads belonging to unique sequences with threshold
I have CLL samples as fastq files and I want to remove those reads which have a unique sequence with less than 10 read counts each. I tried this by following way using awk to make it faster:
for fq in $(ls fastq1/*.fq.gz); do
    basename=$(basename "$fq" | cut -f1 -d '.')
    zcat $fq | awk "/^@/{getline; print}" |sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt; grep -B 1 -A 2 -f tmp.txt $fq | grep -v ^--'| pigz -p 10 > fastq2/$basename.fastq.gz

done

rm tmp.txt

But this is not working well as many reads with more than 10 read counts are also missing. I think it's OK to stick with bash/awk for fast fastq file processing. What am I missing here, or is this command wrong. Can anyone suggest me the right and fast way to do it??

I don't know fastq, but at the start of the pipeline you zcat $fq and in the middle you grep -B 1 -A 2 -f tmp.txt $fq -- do you need to zcat $fq | grep -B1 -A1 -f tmp.txt ?

Either zcat | grep when you're trying to pick reads based on sequences, or zgrep. This code should not work as-is.

@glennjackman: In tmp.txt, I am collecting all unique sequences and finally using this file to count the number of occurrences in the original file. I tried to make it continuous without creating any tmp.txt but that didn't work. Rest other arguments of grep are used to take a complete chunk of 4 lines in fastq files it satisties the criteria of having the number of occurrences more than 10. For more information, you can check the fastq file format.

@RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python.

Your file $fq is either gzipped or not gzipped. If it is gzipped, grep ... $fq should not work. If it's not gzipped, zcat $fq should not work. I don't think both of those statements can work with the same $fq file.

My file $fq is surely gzipped. For this reason, I started with zcat. As per the answer suggested below, zcat|grep can work well with extracted unique sequences in the second step.

You say that you are "trying awk to make it faster". What else have you tried that was slower?

I have not encountered this problem before, and I am wondering if this is an XY problem. Why are you trying to do this?

Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.
I'd adjust your script to the following:
for fq in fastq1/*.fq.gz; do
    basename=$(basename "$fq" .fq.gz)
    zcat "$fq" | awk "/^@/{getline; print}" | sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt
 zcat | grep --no-group-separator -B 1 -A 2 -f tmp.txt "$fq" | pigz -p 10 > "fastq2/$basename.fastq.gz"
done

rm tmp.txt

Instead of

for fq in $(ls fastq1/*.fq.gz)
use this
for fq in fastq1/*.fastq.gz
Iterating over ls is generally not ideal, and can give unexpected results. Better to just rely on bash globbing.

Purely stylistic, but I'd simplify and use basename to remove the suffix, instead of cut.
Rather than removing the separators that grep will introduce by filtering them out, use the following option to not create them:
grep --no-group-separator
And, minor bug, as others have mentioned, you should use zgrep or zcat | grep since your files are compressed.

Oh - and quote your variables.
Finally, while this does work, it may not be the most efficient method. The file is being decompressed twice, and may be better to decompress once to a temporary file that can be used and then deleted.

Thanks for your response. I tried zgrep, but I am not getting any output. Although the script is working till generating tmp.txt. But after that, I didn't see any output file in fastq2. So I tried zcat|grep and it was working. I'll check whether it is giving a correct response or not after completing the pipeline.

@Lot_to_learn  zcat | grep may be a better option depending on your setup, as the --no-group-separator option may not work properly with zgrep

Checking for line start to find the read header is not the best option. Even if unlikely it is possible that the qualitys starts with a @. It's better to compare the current line number. Also you can use awk directly to count and print the result.
zcat $fq|awk 'NR%4==2 {count[$0]++} END { for(sequence in count) { if(count[sequence]>=10) print sequence }}'

As the other statet out, the second command should be fine if you use zcat instead of cat.


text:
How to remove degenerate sequence (N4) at the 5' and 3' of the reads?
I have to analyse GSE114327. For that, I use GSE instruction based on GSM3139597.
In that data processing section has written:

FASTQ reads were subjected to adapter removal (cutadapt) per manufacturer's instructions (NEXTflex Small RNA Sequencing Kit v3).

So, I use cutadapt for trimming. Then

The degenerate sequences (N4) at the 5' and 3' end of the reads were removed

based on the instruction. But I don’t understand that sentence and I don’t know how I can remove degenerate sequence (N4) at the 5' and 3' of the reads.

From the link you posted it looks as though they used UMI tools as per the section which says
'were removed and added to the name of the read as a unique molecular identifier (UMI) using UMI tools'.  The link above has instructions for installation and use.
From this I would guess that this data has Unique Molecular Identifiers (UMIs) for each read.
Looking up the protocol for the nextflex kit used (See Link) confirms this:

Data Analysis
  The 3' and 5' adapters included in this kit both contain 4 random bases that will appear immediately 5' and 3' to the insert in sequencing data. The presence of these random bases should be considered when choosing an alignment strategy. When using “end-to-end” alignment, we
  recommend processing data in the following manner:
  1. Clip the 3' adapter sequence (TGGAATTCTCGGGTGCCAAGG).
  2. Trim the first and last 4 bases from the adapter-clipped reads.
  3. Perform alignments as normal.

It seems that N4 is a 4 base random barcode at both the 5'& 3' end of the read - these are cut from the read and this metadata is added to the read name so that every read can be uniquely identified.
Edit:
You can install UMI tools via:
conda install -c bioconda -c conda-forge umi_tools
or
pip install umi_tools
Usage instructions can be found at:
https://umi-tools.readthedocs.io/en/latest/reference/extract.html

Usage:
  For single ended reads, the following reads from stdin and outputs to stdout:
umi_tools extract --extract-method=string --bc-pattern=[PATTERN] -L extract.log [OPTIONS]
For paired end reads, the following reads end one from stdin and end two from FASTQIN and outputs end one to stdin and end two to FASTQOUT:
umi_tools extract --extract-method=string --bc-pattern=[PATTERN] --bc-pattern2=[PATTERN] --read2-in=[FASTQIN] --read2-out=[FASTQOUT] -L extract.log [OPTIONS]

Dear@Graeme, let me ask that how can I find [PATTERN] for GSM3139597?


text:
semantic similarity measurement for cell line ontologies
I have a set of cell line pairs and I want to know to what extent the pairs are similar based on their ontologies.
The problem I have is that I have found a Python library called Fastsemsim, but it has many similarity measurements. Or to be exact, there are actually many ontology similarity measurements. And I want to know which similarity measurement is best suited for my case (cell line ontologies similarities). There might be a couple of good ones that are good for cells ontologies.
These are a set of similarity measurements that are available in Fastsemsim.
Semantic Similarity

SemSim.TermSemSim
SemSim.ObjSemSim
SemSim.ObjSetSemSim
SemSim.SetSemSim
Specific Semantic Similarity measures
SemSim.MixSemSim
SemSim.avgSemSim
SemSim.maxSemSim
SemSim.BMASemSim
SemSim.SemSimUtils

Source:
https://fastsemsim.readthedocs.io/en/latest/library.html

Welcome to the site. Which is your question. What is the best semantic similarity? Or something else (do you want to compare the cell lines with non-semantic similarities measures)? Do you have other information besides cell lines ontologies (not sure were do you take those from or in which format though)?

The format is OBO but that is not related to the question. My question is which similarity measure (listed) is best or does better at finding semantic similarities between cell lines

About the metric you use, there are several reviews on the topic, if I recall well I think that the information content based are the best ones, of those, I only recognized the JiangConrathSemSim listed on the webpage.
However, if you want to use a semantic similarity measure I would use the BMA method, which is more comprehensive than other methods. But the BMA method is only a way to aggregate several metrics of ontologies into a single numeric value. I am not confident that FastSemSim implements correctly those methods.
Note, that according to the license, you should contact the authors before using it:

You must contact us before using it or its results or any application/tool based on it in any published work.

I am not sure that the cell lines fit in a ontolgogy. Because ontologies are defined as directed acyclic graph (DAG) and there are some cycles that happens in cells like the induced pluripotent stem cells show. However they might be right if you just compare the cell lines and not the state of the cells on those lines (use it only if you want to compare the lineage of your cells lines).


text:
Compare multiple alignment results' aligned bases
I have aligned a nanopore data set to a reference genome with graphmap, minimap2 and BLASR. The alignment results are stored in BAM files.
I would like to do some concordance assessment, looking at the number of base pairs that are mapped to the same position by different aligner.
Are there any tools or ways that allow me to do this effectively?

is this a DNA or RNA dataset you aligned to the reference genome?

it is the RNA dataset from https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md

I don't know of a single tool that can read multiple sets of BAM files and compare them. But you can use the Sam_Reader python3 module to create a csv file that summarizes each position in the genome for each of your 3 results and then compare them manually.
It should be installable with:
pip3 install sam_reader

And then generate the CSV with
from sam_reader import Sam_Reader
my_files = Sam_Reader('my_bam_file_or_folder/')
my_files.per_base_stats(write_file='new_file.csv')

*but FYI this module requires pysam, which does not work on windows

I'm not aware of any tools that do this out-of-the-box. Assuming some programming experience, it would be fairly straightforward to implement this in Python using pysam (or htslib bindings for your favorite scripting language). In particular, the pileup is a very useful data structure.
Ultimately, you'll need to store an index of reference sequence positions by read name and position.
read 1, offset 0 --> chrom 3, offset 248879
read 1, offset 1 --> chrom 3, offset 248880
read 1, offset 2 --> chrom 3, offset 248882
...
...
...

Once you've done this for a pair of BAM files, computing the concordance statistics involves trivial set operations.
There are some important practical considerations though. If you implement the approach completely in memory, it will only be tractable for small-ish BAM files. It wouldn't take very many reads (relatively speaking) to overwhelm available RAM on a laptop or desktop. For large BAM files, you might even overwhelm the available RAM on a big-memory server. Alternatively, you could store the indexes in an on-disk database such as sqlite3. This would keep RAM requirements low, but would slow things down substantially.
I think the most reasonable approach would be to process the BAM files in batches: grab all the reads in a 100kb region (or whatever size makes sense), compute the index, store the congruence statistics, move on to the next 100kb region and repeat, and then aggregate statistics when all regions have been processed.

You can use the bam_alignments_compare.py script from the wub Python package.

Welcome to Bioinformatics.SE, and thanks for your answer. Note that when linking to external sites, please provide some context. Internet links tend to have a short half life, so describing, quoting, or embedding the most relevant content from the link helps preserve the information. See the Provide Context for Links section of this help page for more details.


text:
Hosting IGV web on a server
I have my analysis on a cloud server. The data is quite large, so copying it locally is not really an option. I want to be able to host IGV as a server on that server, so that it can access those files.
I installed IGV web (https://github.com/igvteam/igv-webapp) on my server, and ran it, but the Tracks → Local File option opens a file browser on my local machine, it doesn't let me select files on the server.
Is there any easy way to host IGV web on a cloud server?

Please add some info. Which cloud platform are you using? A commercial solution or an in-house server? What operating system are you using to access the server?

I'm not sure the cloud platform matters. It's just a remote server setup using OpenStack. I'm using Ubuntu both on the server and on my local machine.

It only matters in the in way that you enable X11 forwarding in my solution below. I suppose the specifics of X11 forwarding on different platforms would be better suited for StackOverflow. My answer outlines how to run IGV on a server and access the GUI using an ssh client.

You can use X11 forwarding to run IGV on the headless server where the data is located and simply forward the GUI to your client machine using the X11 protocol. This is a solution I use consistently for my research when working with 20+ tracks that I want to keep on the server.

Download and install IGV for linux from the IGV site on to your server.
Connect to the server using ssh with X11 forwarding enabled
ssh -X user@server
Navigate to the recently downloaded and unzipped IGV folder. Execute the igv.sh script.
./igv.sh

The script launches the IGV session with 4gb of memory. Edit the script to launch with more, 32gb in this edited excerpt.
exec java -showversion --module-path="${prefix}/lib" -Xmx32g \
@"${prefix}/igv.args" \
-Dapple.laf.useScreenMenuBar=true \
-Djava.net.preferIPv4Stack=true \
--module=org.igv/org.broad.igv.ui.Main "$@"

This will probably work the best in my situation. I'm guessing it won't be so helpful for people who don't use X11 on their local machines, though.

My setup also required -Y option, ssh -XY user@server

Also note that this solution need X11 forwarding enabled. Check this question for details.

According to their discription

igv-webapp is a pure-client "genome browser" application based igv.js.

Pure-client, so it will be difficult to let it display tracks that are on the server itself. Perhaps you are looking for something more similar to the Apollo Genome browser. This, while more difficult to install, allows for the display of server-side data tracks and allows for manual curating of annotations as a nice bonus :)

One possibility might be to put the full IGV application on the server, and use the igvR package in R on the server to access the IGV API.  You could then use R plotting facilities etc to show your data, accessed via the IGV API.  There might be APIs for IGV in other languages, if you prefer, not sure.
Edit:  I haven't tried this package, this is a speculative idea.

obviously I'm late but if this is still an open question open an issue in one of our github projects,
github.com/igvteam/igv (Java desktop),
github.com/igvteam/igv.js (pure javascript component),
github.com/igvteam/igv-webapp (html/javascript application, 100% client side)
as appropriate.
As already mentioned the web igvs are 100% client, data is fetched from servers using range requests.   It scales very well,  as long as your data is indexed.   If you really want to generate server side images a better choice might be the UCSC genome browser which is designed for that.


text:
Hard-clip primers in a bam file
I am working on a sample which was sequenced using an amplicon kit and I would like to ensure that no variants are called in the PCR primer sequences. I know that the primer sequences are part of the genome and a variant there is no less valid than a variant anywhere else, but for reasons I cannot control, I need to ensure that no variants can be called in the primer regions.
My first thought was to use a soft-clipping tool like BAMClipper. However, the variant caller I am using (Sentieon's implementation of GATK) preforms a local realignment when calling and if it finds that soft-clipped regions align perfectly against the genome, it will remove the soft clips. This makes perfect sense for cases where the softclips were added by the aligner, but sadly breaks my workflow since I need the softclipped primers to remain clipped.
I would therefore like to hard clip, to completely remove the primer sequences from the bam file before variant calling. I have found one tool that seems to be able to do this, TrimPrimers by Fulcrum Genomics, but this fails on my input bam file with an error I can't really make much of.
Is there any other tool that can take a list of regions (e.g. a bedpe file) and a BAM and hard-clip the regions from the BAM?

Another option would be to adjust the quality of the primer bases to be below the threshold that your variant caller uses.  In my case I have "masked" the primer bases by dropping the quality to 0.  This prevents the variants from being called based that are synthetic but if you have an overlapping amplicon it will call bases from those reads.

That sounds promising. How would I go about that given a bedpe file of primer locations and a bam file? Is there a tool that can take those two files as input and set the MAPQ to 0? Or do you mean changing the quality score in the original fastq file?

The tools we use are custom built but in general we detect the primer for each amplicon and then edit the synthetic (primer) baseQs not the MapQ as 0 for those reads.  You could do it either in the FASTQs prior to alignment or on the BAMs before variant calling.  If you don't have overlapping amplicons this would probably be overkill and filtering the result VCF or redefine the regions of interested to the variant caller might be easier

I know that the primer sequences are part of the genome and a variant
  there is no less valid than a variant anywhere else,

You won't see variants under PCR primers.  You get the primer sequence itself in the vast majority of reads.
I wouldn't clip the bam file, I'd clip either the fastq, or clip within bcl2fastq while making the fastqs.

Thanks, but I do actually get variants in the primer. That's the problem :( And I really don't want to clip the fastq. Both because the primers can help the alignment process and because that's harder to do because I would require fasta sequences of the primers and not regions which is what I have and what assay manufacturers provide. Finally, I never have the bcl files, I will always start from FASTQ (people send me their fastq files).

Variants under the primers can't be real.  Just look at how PCR works.  You will only see a variant under the forward primer if amplification happens from the reverse primer sitting on the orignal template. Once you get amplification happening off of templates made by extending off of the forward primer, the variant is gone.  The orignal template will be outnumbered a billion to 1.

@swbarnes2 primers may be inaccurately synthetized and lead to false positives, theoretically, but I would be more worried about false negatives. A true variant under the primer will be lost, for which you should not call a 0/0 but a ./.

@swbarnes2 that's a very valid point. But my main issue is that I have a client who doesn't want to see any variants in the primers, ever. And I am calling variants in th primers, for whatever reason. It's possible that the primer itself has a mismatch, I guess, and that's what I'm seeing. But I need to completely hide these regions from the variant caller.

Wouldn't it be easier to filter the vcf?

@swbarnes2 indeed it would and that's probably what I'll end up doing. I was hoping there'd be a tool to do this at the bam level to avoid calling them in the first place, but since it looks like there isn't, I'll probably either filter the VCF or extract the target regions from my bedpe file and direct the caller to only call there.

Trimming at the bam stage would altering at multiple fields in the bam; the sequence itself, the quality string, the CIGAR string, the sequence length, plus optional tags at the end, it's a massive headache.


text:
Hierarchical clustering for outlier detection - single cell RNASeq & WGCNA
I wish to conduct WGCNA on a single cell RNASeq dataset and, when choosing the optimal beta parameter after running pickSoftThreshold am presented with a rather strange plot:

I have normalized the data using NormalizeData in Seurat using the LogNormalize method.
One of my colleagues suggested using hclust as a method for outlier detection to see potential cells/samples that could be excluded. I think this makes sense however I haven't found much information on such procedures besides the following paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5845381/.
My question would be, is this procedure just be as simple as running a hierarchical clustering on the data and seeing which samples cluster by themselves (and then removing those based on a minimum cluster size criteria) or is their more to it? Also, for those familiar with WGCNA, why does the scale-free topolgy plot generate such strange looking values, is this indicative of outliers being present in the data?
Thank you for your time.

For outlier identification I suggest using the sample network approach developed by Oldham et al. It basically amounts to constructing a inter-sample connectivity (assuming normalizedData contains the log-normalized data with genes in colunms and samples in rows)
k = colSums(cor(t(normalizedData))),
scaling the connectivities
Z.k = scale(k)
and then removing samples with high negative Z.k (e.g. Z.k < -4).
A slightly different approach is based on euclidean distances rather than correlations,
Z.k = scale(as.matrix(dist(normalizedData)))
and removing samples with high positive Z.k (since distance is a measure of dissimilarity).
Regarding the scale-free topology fit plot, I wouldn't worry about it too much, the fit index sometimes jumps around a bit. It certainly does not look like you have strong outliers or very heterogeneous data. But it is always a good idea to plot the sample clustering tree to check (1) that there are no strong clusters, and (2) that the outlier identification did something reasonable.

Thanks for your response ! I will check out the paper, looks interesting.


text:
Compare functional enrichment between genelists
I want to compare the results from functional enrichment (e.g GO term enrichment generated by gProfileR) between two gene lists. The idea is to see if these gene lists are enriched in the same terms or if they have specific terms.
This is easy to do at a glance when you have few terms but when the number of terms increases and due to redundancies in GO term annotation a manual comparison becomes unreliable.
So the questions is: Is there a tool for comparing the results from functional enrichment between different gene lists?

To compare gene list based on gene ontologies you can use GOSemSim (Note, I'm a contributor to the package). You can select the semantic similarity to compare them, that is how the redundancies and hierarchical relations are counted. The main result you can obtain is similarity matrix for each comparison between gene lists.
However if you perform a GO term enrichment and then you compare the gene list of each GO term by looking at which GO do they belong it won't be much useful. I would suggest you to use BioCor (I'm the developer of it), which can use any pathway information (such as the ones of Reactome) to asses if they are related to the same pathways.

can biocor use output from enrichr or toppgene?

@krushnachChandra I think this is a different question. But yes, you can take a list of pathways from enrichr, toppgeene or other enrichment tools and use it as an input (but you'll need also a source of information about those pathways as input too)

let me try and will update you


text:
For phylogenetic tree construction from core-genome which one is preferable: amino-acid based MSA or nucleotide based MSA?
The genomes are from same species.
Is it true that, in phylogenetic tree constructed from amino-acid based MSA (multiple sequence alignment) some information are lost, so for  phylogenetic reconstruction of closely related sequences nucleotide based MSA  should be used?

For a species complex nucleotide phylogenies all the way.
The reason is neutral mutation, which is not observed at the amino acid level because these mutations are part of the degenerate code which don't induce amino acid substitutions. Such mutations are commonly referred to as the third codon wobble, although occasionally the first codon can mutate and not produce a change in amino acid.
Neutral otherwise known as synonymous mutations occur at a much higher rate than mutations that cause a substitution of an amino acid. This is the reason why  you were told information is lost at the amino acid level because you're losing neutral mutations.
There are circumstances when the amino acid phylogenies are preferred to nucleoside trees, but these are bit pretty specialised areas.

also nucleotides include the (often variable) noncoding sequences that don't encode aa


text:
Paired end sequencing with R1 and R2 of different length. Possible?
There are commercial sequencing kits/sequencers that allow for paired end sequencing in which the two reads obtained for each fragment are of different length?
For example 150bp for R1 and 50bp for R2?
Is it at least possible in principle? There are some published speculations about that possibility?

it's used for certain applications, e.g. 10X single cell where R1 is shorter, being only the bar code and UMI

If you're looking for technologies that do this explicitly to optimize for something, it's possible in principle. If you're looking for an explanation why some files you're looking at are different, it may be due to some pre-processing on the files, such as in the scRNA-seq case

Yes. For Illumina kits for example you include R1 and R2 length in the sample sheet. As long as you don't exceed the maximum number of cycles you will be fine. For a 150 cycle kit you can perfectly run 75-75, 125-25 and even 50-10.

As indicated in the comments, yes, you can definitely tell standard Illumina sequencers to sequence mates in a pair to different lengths. This is quite common in single-cell RNA-seq, where R1 is often just a cell barcode and an UMI, so sequencing it to 150 bases would be a huge waste or reagent. One caveat to this is that you can't buy reagents sized for these custom length, just things like 150 base PE or other standards. So unless you end up saving and pooling reagents after such runs you end up paying the same for projects with one shorter read than if you'd just directly sequenced everything to the same length. Given that sequencer runs only take a day or two at the moment, this also doesn't result in that immense of a time savings.


text:
In SBOL, when should I annotate a DNA sequence vs. making sub-components?
When representing a DNA sequence in SBOL as a ComponentDefinition, you can mark things like promoters and coding sequences in two different ways:

either as a SequenceAnnotation in the ComponentDefinition, or
as a Component linking to a ComponentDefinition for the sub-comoponent

When should you choose one over the other?

@James McLaughlin. If there is some error on the question and answers about the class I think it would be better if you provide a comment it could be a misunderstanding (unless it has changed and you are updating the question which it would be nice to know in the edit summary)

@jakebeal, there is a pending edit for both your question and your answer (you might not have access to those links, if you don't just click on the "edit" under each or the notification in your inbox) changing SequenceFeature to  SequenceAnnotation. Please confirm which of those is the correct class name.

@terdon Ah, James was right --- I typo'ed the names.  Edit accepted.

Typically you want to use a ComponentDefinition when an sub-structure is something you might want to pull out and re-use in another genetic design, and a SequenceAnnotation whenever it is not.
For example:

Promoters, terminators, and protein coding sequences are typically best represented with their own ComponentDefinition, since these are often pulled out and re-combined in new genetic designs.  Representing such an element with a ComponentDefinition makes it easy to find the commonalities between designs.
Assembly scars, binding sites, and cut sites are typically best represented with a SequenceAnnotation because they are typically not very interesting except in the context of their surrounding genetic context. Representing such an element with a SequenceAnnotation helps keep them associated tightly with that context.

Not every case is cut and dried, however.  For example, ribosome entry sites are often treated as separable components, but their performance is very tightly tied to the specifics of the coding sequence that they modulate.  Thus, in some cases it may make more sense to represent them separately as a ComponentDefinition and in others to represent them as a SequenceAnnotation on a joint RES/CDS ComponentDefinition.

Note that this answer was written for SBOL 2.  For SBOL 3, change ComponentDefinition to Component and SequenceAnnotation to SequenceFeature. There are also shortcuts in SBOL 3 for representing things in other databases like UniProt and ChEBI


text:
Mapping GPL Platform Number to Microarray Name
I am trying to download a dataset from GEO using the GEOquery R package. Using this package, I can retrieve the GPL platform number(s) for the dataset, and I hope to use a CDF file from BrainArray on the dataset in order to undo any existing normalization performed, so that I can have multiple datasets that are normalized in a standard manner.
In order to do this, I need to know what the microarray platform name is. This information is usually provided on the website online (for example, here, you can see that next to GPL96, it says HG-U133A, which is the CDF I would use from BrainArray).
However, I am wondering if there is any way to parse this platform name using the GEOquery package in R itself rather than having to manually check for every dataset and compare it to the BrainArray website?
If it is not possible to do this in the GEOquery package, is there any online tool for mapping GPL to platform name? I found this online which has the mapping for a decent amount of GPL numbers, but not all of them, meaning it sometimes crashes. I was wondering if there is a foolproof way to do what I am trying to do.
Thanks for the help in advance!
Edit:
I know I can access the highlighted portion in the first image below (the GPL platform number) through GEOquery. I'm wondering if I can access the name of the platform (highlighted in the second image). If not, my alternative option will probably be to scrape the web page to find this information, but if its possible to do so through GEOquery I'd rather do that.


text:
Does SBOL require representing intermediate products like mRNA?
Do I need to include intermediate gene products like mRNA in order to express a gene regulatory relationship in SBOL?
For example, let's say that I'm representing transcriptional regulation of the pTet promoter by means of a coding sequence that expresses TetR. Do I need to represent the TetR transcript mRNA as well as the TetR protein?

As far as I understand, the purpose of SBOL is not to create a comprehensive model, but to communicate a design.  So the question you have to ask is "is the mRNA part of my design, or is it an implicit part of the biology which does not help to convey what I'm trying to build?".
I think your summary is a giveaway:

For example, let's say that I'm representing transcriptional regulation of the pTet promoter by means of a coding sequence that expresses TetR.

Notice that the mRNA was not mentioned?  I also daresay that if you drew this design as a diagram, the mRNA would be omitted - another indication that it is not helpful in communicating what you designed.


text:
Why does SBOL have both Component and ComponentDefinition?
What is the point of having both a Component and a ComponentDefinition in SBOL? Shouldn't you be able to link straight from one component to its sub-components?
The same question goes for Module and ModuleDefinition.

The key to understanding this is to think about a system with more than one copy of a genetic construct.
For example, let's say you have a construct that uses two copies of the pBAD promoter, each regulating a different coding sequence. Both of these promoters are pBAD promoters, with the same sequence and type information, so they should share the same ComponentDefinition.
But there are two copies. In the context of the larger construct, we need to be able to talk about each copy separately. For example, we might need to say that the first comes after a MoClo cut site, while the second comes after a terminator. For that, we need the Component object, to have a way to say "this particular copy of pBAD", as opposed to "pBAD in general", which is what the ComponentDefinition tells us about.

Essentially, a ComponentDefinition represents the component itself, while a Component represents the component being used as a sub-component.
There is a proposal within the SBOL community as to whether for SBOL version 3 ComponentDefinition should be renamed to Component and Component should be renamed to SubComponent to make this situation less confusing.


text:
How do I select a proper design matrix to use in DRIMSeq?
I am trying to use DRIMseq for DTU with 2 treatments on two different strains of an animal model.
How do I select a proper design matrix to use in DRIMSeq ? How I know that one is the correct one to use? And how do I see which transcripts are differentially used among conditions after the pipeline was run? I suspect it is in the plotProportions function.
Also, where do I see the logFoldChange between conditions?
Many thanks in advance as it would help me a lot in my work.
> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sva_3.30.1          BiocParallel_1.16.6 genefilter_1.64.0   mgcv_1.8-28
[5] nlme_3.1-137        edgeR_3.24.3        limma_3.38.3

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1           compiler_3.5.3       pillar_1.4.2         bitops_1.0-6
 [5] tools_3.5.3          digest_0.6.20        zeallot_0.1.0        bit_1.1-14
 [9] annotate_1.60.1      RSQLite_2.1.1        memoise_1.1.0        tibble_2.1.3
[13] lattice_0.20-38      pkgconfig_2.0.2      rlang_0.4.0          Matrix_1.2-17
[17] DBI_1.0.0            rstudioapi_0.10      yaml_2.2.0           parallel_3.5.3
[21] IRanges_2.16.0       S4Vectors_0.20.1     vctrs_0.2.0          locfit_1.5-9.1
[25] stats4_3.5.3         bit64_0.9-7          grid_3.5.3           Biobase_2.42.0
[29] AnnotationDbi_1.44.0 survival_2.44-1.1    XML_3.98-1.20        blob_1.2.0
[33] matrixStats_0.54.0   backports_1.1.4      splines_3.5.3        BiocGenerics_0.28.0
[37] xtable_1.8-4         RCurl_1.95-4.12      crayon_1.3.4

Welcome to the site. The design depends on the question and on the conditions of the experiments. It would be easier for us to help you if you provide a table with the information about each sample sequenced. The strain, treatment,  sex,  batch and other relevant conditions used of each sample. Also it might (or not) be relevant to know the tissue you are studying (or if you are mixing several tissues on the analysis)

From your question I understand that you have two treatments in two mouse strains.
If you are interested in general treatment effects independet of the animal strain you would use the both variables in the design formula (similar to controlling for a batch effect).
meta.data <- data.frame(treatment = c("a", "a", "b", "b"), strain = c("x", "y", "x", "y"))
design_combined <- model.matrix(~ treatment + strain, data = meta.data)

If you want to make the comparisons for each strain separately, the easiest way will be to split your data set by strain and do a set up just for treatment. Alternatively, you can also merge the two columns into a new one in case you also want to also compare between the strains.
meta.data$combined <- paste0(meta.data$treatment, "_", meta.data$strain)
# refactor to change the reference level
meta.data$combined <- factor(meta.data$combined, levels = "bx", "by", "ax", "ay")
design_combined <- model.matrix(~ combined, data = meta.data)

To get the fold changes, you can extract the proportions for each transcript using
DRIMSeq::proportions(dataSet)

And from the proportions the log2 fold changes can be calculated with ease.
However, when I look at your plot something seems very off with your analysis right now. In all the examples I have seen and also in my own data the diamonds are representing the mean proportion of a transcript in a condition. See this example from the DRIMSeq manual on Bioconductor:

In your case there are several different means per group and if you look at the line to the very right the two treatments seem to be to be right on top of each other.
So I think that treatment and strain are somehow mixed in your setup and this is causing this kind of merged plot. If you use a design as outlined above you should get much better results. Additionally, you could maybe check that the meta data is correct.
Also in case you have not read it yet, the manual is a great resource:
https://bioconductor.org/packages/release/bioc/vignettes/DRIMSeq/inst/doc/DRIMSeq.pdf

How do i calculate the LogFC from the proportions? It isn't in the manual.

divide the proportions and take the log2. Log2(proportion A / porportion B) = Log2FC of the proportions with the condition B as the reference level.

Thank you very much for your responses, they are helping me a lot. When i do the command  DRIMSeq::proportions(dataSet) , I get proportions for each transcript for all samples. Should I sum the proportions in the samples that belong to category A and the same for samples that belong to category B, and then finally do Log2(proportion A / porportion B) ?

Glad to help. Yes, you can calculate the overall proportion from the single proportions but you should calculate the mean proportion for each experimental group (which is the same as the diamonds in the proportion plot) and then calculate log2(mean(propA)/mean(propB)). Otherwise you run into problems when your sample sizes are not equal.

Very clear! Thank you. I'm now computing the DRIMSeq pipeline for every possible combination of the design matrix I could think of (they are 8) and see any give results like the ones you described above. Despite being computacionally expensive do you think its a good idea, in order to select a proper design matrix ?


text:
making a customer reference by using cell ranger
Anyone here knows how to add a new gene (GFP) to customer reference by using cellranger? We use the Cre/loxp system by using R26RmTmG mice and want to use GFP to isolate specific cells. I have read the cellranger workflow but still, have no idea how to do it. Like it is said that I should update the GTF file as shown in the screenshot.

Step 1:
Add the the GFP sequence to your fasta file
>GFP_name
ATG.....TAA

Step 2:
Annotate the gtf so that it contains all the necessary information to annotate your information
Add this line to your GTF file but modify where necessary:
"Make sure that all columns are tab separated"
GFP_name    AddedGene   exon    1   1966    .   +   1   gene_id "GFP"; transcript_id "GFP"; gene_name "GFP"; transcript_name "GFP"

The first column has to correspond to the header of the fasta sequence
start and end indicates the length of your sequence. I added a 1966 bp long sequence to my gtf

Code example:
step 1:
cat reference.fa gfp_seq.fa > new_reference.fa

step 2:
cat genes.gtf gfp.gtf > new_genes.gtf

-  gfp_seq.fa is your custom gene in fasta file
-  gfp.gtf is your custom gtf line or lines if you want to add more
You can also open your current reference fasta and gtf file in a text editor and copy paste the sequence and gtf line in there.
Then run cellranger mkref on your new_reference and new_genes file

Hi Mack, thank you for your reply. And I was still wandering whether "the genes.gtf" in Code2 "cat genes.gtf gfp.gtf > new_genes.gtf "  is the gtf that had to be filtered before or after Code 2 to proceed filter step?

You can add this line to the filtered gtf file. You can download the gtf file for your species of interest, format it for cellranger using the mkgtf function. The output.gtf will be used together with your reference.fa. You can add all custom genes to this output.gtf as long as you follow the correct format and annotate it as exon. You can also add gtf lines to the non-filtered gtf file and filter it then. I prefer to generate a species reference for myself and the lab, and modify that one for my own use.

Hi Mack, much appreciated for your help. And after making reference, it comes the error as below: " Property 'transcript_id' not found in GTF line 1704071: GFP AddedGene exon 1 846 . + 1 gene_id “GFP"; transcript_id “GFP"; gene_name “GFP"; transcript_name “GFP" "

I am sorry for that. I believe I had a copy/paste error that adjusted the quotes. Make sure your quotes around the "GFP" are correct in your gtf line. That should solve the problem but feel free to reach out if not. Also make sure each column is tab separated.

Hi Mack, thank you for your advice. Finally, I make it. And I have another question want to ask you whether it is possible that I can add a truncate sequence into the customer reference? The information is as below: 1. We use Cre/loxp system. 2. The truncate we want to add into customer reference is Nkx2-5.  Which was excised of exon2 and the flanking sequences resulted in a truncated RNA molecule that lacks exon 2 coding sequences, the splice sites, and the 3' noncoding region, rendering it essentially nonfunctional.

And maybe you can have a look at its construction in the paper. Here is the link of that mice https://ars.els-cdn.com/content/image/1-s2.0-S0092867404004052-gr1_lrg.jpg. You can see the Fig.1A.

Technically yes, but if the transcript is stable and detectable by the 10x chemistry, you should be able to map it to the wt-nkx annotation. You might need transcript assembly instead of mapping to a reference but that is material for a new question.

Thank you for your time. And the last question is whether we add the GFP into the reference will make the result slightly different?

That depends on what you are trying to assess. If you post an analysis with and without GFP in reference and explain the question you are addressing it will be more helpful for us to help you out.

It shows that the web_summary information changes a little bit. Like before using GFP reference, the cell number is 6542, the mean reads per cell are 35877, and the median genes per cell is 2539 all of them changed to 8452 cell, 27770 reads per cell, 1559 genes per cell. So did that cause any problem? Did that usually happen or I did something wrong?

Even tho cellranger is intuitive there is still a lot of bioinformatics that happens under the hood. Without knowing how you ran cellranger it is hard to say if this is supposed to happen. For example: are you sure you ran both references on the same cellranger version? Version 2 and 3 use different algorithms to call cells. Did you enter the same parameters in cellranger count (--expect-cells vs --force-cells vs default, etc?). A new question with details on how you ran cellranger would be helpful to troubleshoot or explain what happend.

Hi Mack, thank you for your advice, after looking my cell ranger version it is 2.1.0 before making GFP reference. And I change it to 3.0.1 they work right now!

Update your gtf file, update the genome fasta, run cellranger mkref.

This looks like a comment. Or at least I don't understand how this helps to add a new gene. Could you expand it so that it is more clear ?

Although this might be a valid answer, I find somehow dissatisfactory. I think some explanation would improve it a lot.


text:
Extract mapping coverage from GTF files
I am doing paired end RNA-seq analysis. I used STAR mapper and then StringTie for getting quantification of my data. After successfully running StringTie, I got .gtf files which will give me information about genes, their coverage along with TPM and FPKM values. It also produced some more .ctab files.
I want to know how what percentage of my reads mapped to reference genome. I know this can be done with BAM or BED files produced after STAR alignment. However, I want to know if this can be done with output files generated by StringTie?
I just want to clarify my question after reading comments and answers I am getting. I want to know if this can be done with .gtf files of StringTie. As I mentioned before, I already know how to achieve this with BAM files. This is little bit weird situation because, I have deleted my BAM files to get more space on my hard drive. So I could just run entire analysis again, but wanted to know if it is possible without redoing whole analysis. Apologies for the confusion.

So consider this, you know how many reads map to each of your genes (1). Does that tell you anything about (2) the number of reads that mapped against non-gene regions (3) the number of reads that did not map against the genome at all?
If you only use the gtf files, you know the number of (1) but not of (2) and (3). To me it sounds like you would want to do something as described in this post.
What it boils down to is:
Finding out how many (unique) reads are mapped, counting pairs once
samtools view -F 0x4 foo.sorted.bam | cut -f 1 | sort | uniq | wc -l

Counting how many total reads you have. Since you have paired-end data the number of forward reads suffice:
cat forward_reads.fastq | grep '@' | wc -l

The percentage is then simply mapped_reads/total_reads

@terdon sorry about that, I updated my answer

Perfect! Thank you and welcome aboard!

I want to know how what percentage of my reads mapped to reference
  genome.

And the Log.final.out file from STAR doesn't tell you that?


text:
How do you represent small molecules in SBOL?
In SBOL, how do you represent small molecules like arabinose or chloramphenicol?
There's a type for saying a component is a small molecule, but they don't really have sequences like DNA or proteins do, so what's the recommended way to provide information about it other than just the name?


text:
Single Cell RNA seq Analysis for low input Cells?
I am having a single cell RNA seq data from ~500 Cells. I do understand this number is very low for present day analysis tools like Seurat. So, I want to know what if one have a such a small number of cells which tool works best?
Thanks

you can perfectly use seurat as part of your analysis pipeline. It is hard to suggest tools or statistics without knowing the type of data you are looking at and what questions you want to address.

@Mack123456 Thanks for your comment. I am trying to use single cell RNA seq data on MiNION. On a trial basis I just started with ~500 Cells.

Hi, you will need to forward as much information as possible, or simply accept @Mark123456 answer. The problem is that we don't know the objective of the study.

@MichaelG. Thanks. I am ready to accept the comment as answer. But I cant find the accept answer tag in comment section.

Hi @user3377241 no need these are just comments not a full answer. If you look at the other questions it's detailed technical advice rather than broad answers. If you have a specific technical issue you can return and ask about it

You can perfectly use any of the single cell processing packages out there. Seurat, Monocle3, Scanpy, Kallisto, etc. Following some of the tutorials and vignettes will most likely help you figuring out how to proceed.
These packages provide useful functions to explore single cell data, process the data and each of them has a series of advanced features to address specific questions. You don't need 10k cells to look at number of gene/cell, transcript/cells, % mitochondrial genes, etc. You can perfectly normalize and standardize 500 cells and run PCA, and analyze the complexity that describes your data.
To address if 500 cells is enough for you pilot we need to know what question you want to address.
The following tool for example gives some intuition on how many cells to sample if you wish to detect n-cells per cluster.
https://satijalab.org/howmanycells
Hopes this helps