text:
HMM Profile from Convergence of retrotransposons in oomycetes and plants
I am trying to recreate the experiment from Convergence of retrotransposons in
oomycetes and plants by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. In the "Methods" section a

HMM profile was constructed from the amino acid alignment of aRNH
  sequences from the Ustyantsev et al

I was wondering where I could get this profile. I checked the cited papers but they were also referencing other papers. Could someone maybe clarify where I could get the profile?

Your question is how to to this HMM or how to obtain/download it? In the later case you can always try to contact the authors if it is important to use exactly the same HMM profile

The alignment they used to build the HMM is included as supplementary file 8 in the paper you cite:

Multiple amino acid sequence alignment of RNH genes and domains from diverse taxa constructed and used for the phylogenetic reconstruction in the present study. (TXT 45 kb)

So, you can download that alignment and then use a tool like hmmbuild (part of HMMER) to build the HMM profile:
hmmbuild model.hmm 13100_2017_87_MOESM8_ESM.txt

The command above will create the file model.hmm which contains the HMM built from the alignment in the supplementary file 13100_2017_87_MOESM8_ESM.txt. I am including the first few lines here (this is not the entire file, don't use this as the model—build it yourself. I only include it to illustrate the format):
HMMER3/f [3.1b2 | February 2015]
NAME  13100_2017_87_MOESM8_ESM
LENG  121
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Mon Dec 11 13:33:29 2017
NSEQ  165
EFFN  5.697556
CKSUM 2893352190
STATS LOCAL MSV       -9.8975  0.71332
STATS LOCAL VITERBI  -10.6450  0.71332
STATS LOCAL FORWARD   -3.9247  0.71332
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.34473  4.21240  2.85435  2.69203  3.54228  2.79153  3.54812  3.01140  2.72630  2.48370  3.82225  2.94269  3.49156  3.05882  2.84632  2.63343  2.86306  2.72633  4.27012  3.56181
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.13773  5.45955  2.08419  0.61958  0.77255  0.00000        *
      1   2.99488  4.57852  4.44659  3.29307  1.80597  4.16203  3.88788  2.87265  3.75733  2.38594  3.30344  3.33633  4.53284  3.64324  3.69630  2.90350  3.30742  2.83413  4.36991  1.50673      1 y - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00723  5.32905  6.05140  0.61958  0.77255  0.45275  1.01027
      2   2.66602  2.17750  5.01715  4.40788  1.72617  4.29647  4.62737  2.89579  4.21116  2.71668  3.64495  4.43195  4.66188  4.35829  3.99238  3.10363  2.02207  1.89528  5.14692  3.02155      2 f - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00661  5.41881  6.14115  0.61958  0.77255  0.61842  0.77391

thanks for the clarification. I missed the additional file since I was looking at a printed out copy.


text:
How to use the hmmsearch on prototypic repetitive sequences Repbase Update database?
I want to use the hmmsearch proposed in Convergence of retrotransposons in
oomycetes and plants by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. After the help from terdon I managed to create the correct profile I assumed that the command
hmmsearch aRNH_profile.hmm prirep.ref > output.out

would be the correct command. I downloaded the compressed database RepBase22.11.fasta.tar.gz from here in the fasta format.
However when unpacking I get the following tree
.
├── angrep.ref
├── appendix
│   ├── angapp.ref
│   ├── athapp.ref
│   ├── bctapp.ref
│   ├── cbrapp.ref
│  ....
├── athrep.ref
├── bctrep.ref
├── cbrrep.ref
...

which file is the correct fasta file? I assumed prirep.ref  might be correct file but the output shows
[No hits detected that satisfy reporting thresholds]

The site requires an account in order to download. Isn't there a README or similar file explaining what each file contains? What other files are there?
Are they all .ref?

true. I created a account. All the files are .ref files and there is no README file

OK, so what species are you looking for? All plants? Plants and fungi? Only fungi? How about algae? What is the "correct" file for you?

RepBase explain what these files are on their website. It is usually a good idea to check the documentation when downloading files so you know what you're downloading.
Based on the list of files described on the linked page, there is no prirep.ref file, so I assume you mean plnrep.ref which is "Other Plants". In any case, again according to that page, the files with plant and/or fungal data (which is what I assume you want; you never said) are:
athrep.ref  Arabidopsis thaliana
oryrep.ref  Oryza sativa
dcotrep.ref Other Dicot Plants
mcotrep.ref Other Monocot Plants
grasrep.ref Other Grass Plants
plnrep.ref  Other Plants
fngrep.ref  Fungi
There's also chlrep.ref for the green alga Chlamydomonas reinhardtii, and diarep.ref for the diatom Thalassiosira pseudonana. I don't know whether you would want to consider either here.
In any case, since these are fasta files, you cna combine them by simply cating them together:
cat athrep.ref oryrep.ref dcotrep.ref mcotrep.ref grasrep.ref \
    plnrep.ref fngrep.ref > plantsAndFungi.ref

The, run your search on the plantsAndFungi.ref file.

as I was talking to the author he mentioned that they concatenated all the files and ran a translator on nucleotide on the data. Before running the search

source: Stack Exchange

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

sort: votes


text:
LD analysis in PLINK based on reference and a SNP list
I'm new to PLINK and genetics, and getting confused with two PLINK commands for LD analysis:

plink --bfile hapmap --r2 --ld-window-r2 --ld-snp-list --ld-window
plink --noweb --bfile --clump hapmap --clump inputfile --clump-field P --clump-p1 --clump-p2 --clump-r2 --clump-kb

It seems to me that for #1, the input files are a reference data and a SNP list file, while for #2, it is a reference data and a .assoc file.
In the output for #2, p-values are included with other results (according to PLINK manual) which is not the case with the output from #1. If I'm interested in finding proxy SNPs for the SNPS on my list in hapmap data, will it be sufficient for me to use the output from #1?
After invoking command #1 using PLINK 1.9, its output gave me 6 columns: "Chromosome number", "chromosomal position", "SNP name" for query/target SNP and 1 column of "R2". I cannot carry out command #1 because I do not have the .assoc file and only have a list of SNPs.
My main question is - what is the difference between these two commands and which of the two should give me the best result for finding proxy SNPs for my list of SNPs.

It looks like the first command is calculating pairwise linkage disequilibrium r2 between your SNPs.
The second command is clumping your dataset based on a p-value over a file of association test results, which I guess is a pre-cursor to computing polygenic risk scores.
Clumping is the process of thinning your SNP list to only keep uncorrelated SNPs so that each representative region has only one marker. In this case you are asking PLINK to keep the marker with the highest p-value.
I'm not sure what your definition of proxy SNPs is but if you are looking for tag SNPs to a SNP of your interest, you should go with command 1 and choose a proper LD r2 threshold.

Many thanks @Vivek for your explanation. I really appreciate it!

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


text:
low-memory high-speed replacement for Picard MarkDuplicates
I am running Picard MarkDuplicates with the following parameters below. On the file described, it takes about 41.6Gb of RAM memory and about 20-25 minutes to compute (only uses 1 core AFAICS).
java -Xmx54890m -jar picard.jar MarkDuplicates \
    I=/home/dnanexus/in/sorted_bam/TST73-86-2IC_S16_L00.bml.GRCh38.karyo.bam \
    O=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.deduplicated.bam \
    M=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.duplication_metrics \
    CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT REMOVE_DUPLICATES=true

What equivalent free/open source software (samtools?, sambamba?) can I use to obtain the same deduplicated bam file but at a lower cost, considering the instance used here has 8 cores but only one of them is being used, and has about 68Gb of RAM out of which it's using 46Gb odd at peak.
I think previous versions of samtools/sambamba didn't behave the same way as Picard MarkDuplicates with the parameters above, but I would like to know which recent versions would give me the same output.

I assume you only want free (in terms of price, at least) software, right? Sentieon have a tool for this that is orders of magnitude faster (and I would guess uses less RAM but I'd have to check) but you need to pay for it.

Correct, open source.

I done no benchmark, but I guess samblaster is bit faster (by my subjective perception of time).

You should specify a smaller "-Xmx". MarkDuplicate doesn't need that much memory, but seeing a large "-Xmx", java will use extra.

As you suggested, sambamba is faster at marking duplicates than picard (it's also multithreaded). Recent versions of samtools have a rewritten duplicate marking algorithm, though I doubt it'll be as quick as sambamba. 46GB of RAM seems excessive for marking duplicates unless you're having it store the whole file in memory.

Does sambamba with default parameters equate to Picard MarkDuplicates?

I'm not sure that any tool exactly equates (there are a lot of edge cases in duplicate marking and picard misses some things), but it's a sufficient replacement. Sambamba is also used for this purpose by a number of international consortia (I think either DEEP or IHEC, if I remember correctly).


text:
AND or OR gene - protein logical relation in KEGG
How do I know if it is a AND or OR relationship between genes associated with the same enzyme in KEGG? E.g., for EC:1.6.1.2 (http://www.genome.jp/dbget-bin/www_bget?ec:1.6.1.2) does the bacterium Streptomyces coelicolor (SCO) need SCO7622 AND SCO7623 or SCO7622 OR SCO7623?

The enzyme EC:1.6.1.2 is an NAD(P) transhydrogenase and consists of two subunits. The two subunits are encoded by two different genes pntA and pntB.
So in terms of computer logic I would say AND, both subunits are necessary for a functional NAD(P) transhydrogenase enzyme (EC:1.6.1.2).

This was only an example, is there a way of knowing in general?

With enzymes from KEGG there is not, I am afraid. In your example it becomes clear that we see 2 subunits, and hence form the same protein. However, if you look for example at human isocitrate dehydrogenase (IDH) EC 1.1.1.42, you will see that two different genes (IDH1 and IDH2) are listed. These 2 genes are not subunits of the same protein but are 2 different proteins (one in mitochondrium and the other in cytosol). With KEGG the interest in enzyme, not gene.


text:
Protein fold pathfinding?
This is probably a long-shot since it would be so dependent upon the underlying folding algorithms. Are there known algorithms for mapping the path followed in a protein folding simulation when you have the starting conformation and the solution conformation such that you can count the number of steps involved in going from one to the other?
The starting conformation would be something generated from MSM Builder while the solution conformation would be the result of hours, days, weeks, months, etc of folding.

Wouldn't all algorithms do that? I mean, if the algorithm is following a path, you would just have to implement it so that it prints out a number for each path. I don't see how this could be algorithm-dependent. Did you mean to ask whether there are any implementations, any tools implementing an algorithm, that already do this?

One way to get one of the many statistically possible paths is to start from the end and work backwards. This works with the maze puzzles on kids' menus, but also with proteins. To find how a ligand enters an active site, dynamic undocking uses tethered MD to pull a bound form out. Likewise for a protein folding starting from a folded state and lowering the attraction term of Lenard-Jones or gradually increasing the reference temperature used by the Berendsen thermostat in an MD run (as done for DNA annealing simulations).
This assumes that the system is reversible (not true for protein with a maturation step like GFP), but it is avoids getting stuck in local minima or dead-end paths in the case of maze puzzles.
It also assumes that peptides start folding only once they are finished being synthesised, which is not true and actually circular permutations are an interesting case as many the permutation windows/spokes are unfolded.
However, the path down a folding funnel is more of a hike/roll down a canyon than a maze puzzles on paper as there are many possible routes to the solution. So having a bunch of trajectories makes a nice video, but does not say much...

A seemingly similar concept, but very different, is the pathway reconstruction for the evolution trajectory that leads to a multi-mutation variant of interest. This is well studied because the intermediates are not transient but have to persevere in the population until another mutation appears. If they have a too high fitness cost they will not be around long enough for a compensation mutation (cf. epistasis).

Answer from @Llrs, converted from comment:
I would say that the problem is defining each steps on your folding path, it is a continum, I don't know if you have seen one of those videos/gifs of a protein folding, they are quite illustrative on this aspect. Usually there isn't a clear step with a pause unless it reaches a local minimum, it could be possible to save that minimums with its folding but usually they take time and resources to calculates the folding process for a long time, so probably one ends up with just a few steps of the folding.

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


text:
CRISPR Sequence Finder and Database Download
I am searching for tools to pull CRISPR Spacers from Bacterial Genomes. I am aware of the CRISPRDB and the corresponding identification tool on the web server.

Are there other tools for finding CRISPR Spacers?
Is there a way to download the spacer sequences from CRISPRDB?
Is there another database that contains all validated / inferred bacterial CRISPR sequences?

CRT (CRISPR Recognition Tool) http://www.room220.com/crt/
Crisprs Finder Online Tool https://crispr.i2bc.paris-saclay.fr/Server/
Piler-CR http://www.drive5.com/pilercr/
(Omic Tools List) https://omictools.com/crispr-detection-category


text:
Is it possible to check if a patient has the HLA-B27 antigen from his exome stored in a VCF file
Is it possible to check if patient has the HLAB27 antigen by comparing his exome stored in VCF with all known HLAB27 variants?
For instance, rs4349859 is a known HLAB27 SNP.
If I can find rs4349859 in the patient's VCF file, could I say he is HLAB27+ ?

What do you mean "test" it? VCF files normally don't have nucleotide sequences (well, they have the variant, but that isn't what you normally need). If you want to check whether a patient has the antigen, you should probably be looking at protein data, not DNA. Please [edit] your question and clarify what you are trying to do.

Some papers claim that you can do HLA typing from exome sequencing data. I have not validated this myself. I also guess it depends on capture. Nonetheless, I don't think you can type HLA with VCF generated with a standard pipeline, which is not specifically designed for HLA. Current variant callers won't work on HLA genes due to their high diversity. Also in general, there are thousands of HLA-B alleles. In most cases, you can't find the allele by looking at only one SNP.
No, you can't check if the patient has HLA-B27 only from VCF. There might be a hope if you get raw reads and run a proper HLA typing tool designed for such data.


text:
How to rebuild the workflow for mining in oomynete genomes in UGENE
I want to recreate the UGENE workflow proposed in Convergence of retrotransposons in oomycetes and plants. After being able to understand the data and hmmsearch Here and Here and getting most of the result I want to try to recreate the workflow in UGENE.
I manged to implement the workflow for the search
However in the latter parts I struggle to implement the steps of the algorithm

sequences surrounding the regions of significant similarity to the aRNH profile were expanded, when possible, to 10,000 bp in both directions

How can sequences expanded in UGENE?
How can implement the next steps also? for step three

Third, the enlarged sequences were screened for the presence of significant similarity to RT domains of non-LTR-RTs and LTR-RTs HMM profiles using hmmsearch.

I assume it is very simmilar to the first step in building the right profile and running hmmsearch but the following steps are still unclear to me.

Expanding means looking outside the original sequence (in this case by 10kb), so you need to find the position of the sequences and look in the genome for a larger range. Could you clarify what problems do you have with the step three? To me it means to use hmmsearch with specific profiles to search in those sequences

Well I thought that the search result of step one returns the loction of aRNH  then from that location expand in the old data by 10k bp and hmmsearh with the new LTR-RT and non-LTR-RT profiles. but how is the exapnsion from the search result done that was unclear to me.


text:
How to get all PDB homologs from Uniprot (mapping + BLAST)?
I'd like to create a dataset consisting of all sequences which are either present in the PDB, or whose homolog is present in the PDB. In other words, any sequence in the PDB or any sequence related to it. The similiarity margin is to be very wide, so anything above, for instance, 25% sequence identity is accepted (perhaps there are better criteria, but you get the point).
No further limitations, so any organism, with or without experimental data, not limited to the manually curated Swissprot set.
I have tried the SIFTS database, but that contains, as far as I can tell, only mappings for sequence identity >85%. I would like to go far below that.

Thanks for the edit, that's much clearer! So, next question: how do you define "homologs"? Are you looking for functional homology alone? Is sequence homology enough? I am guessing you will be assuming that homologous sequences will share a structure but that is most certainly not always the case. Is there any reason you don't do it the other way around? Get all sequences from the PDB and then map those to UniProt accessions?

@terdon Homologs are defined here as >25% sequence identity, so very very broadly. I would, in fact, be willing to run a BLAST search for each PDB entry vs Uniprot, but I assume this will take quite some time, so I'm checking to avoid reinventing the wheel.

For the correspondence between PDB and Uniprot entries you can use SIFTS -- a semi-automated mapping between PDB and UniProt maintained by PDBe.
The pipeline that creates the mappings uses BLAST and a few other criteria to decide which UniProt entry should be assigned to each PDB entry. The SIFTS website has all the data in CSV files.

Thanks for the advice! I took a look. In the referenced publication, I only found a mention that they are looking for a 90% sequence identity to accept a cross-reference; this is not enough for me - I need all UniProt entries with, say, sequence identity above 25%.
As I understand it, this is what is contained in the files on the ftp server ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/ , notably the uniprot_segments_observed.

You can download a BLAST db containing all sequences of proteins in Uniprot and in the PDB. The way I would go about this is first download the databases for uniprot and PDB, then query the PDB database for each sequence from Uniprot. If you get a BLAST hit above a certain threshold (whatever you define as homolog) then add that sequence to a file, and voila, you have a full set of sequences that fulfill your criteria. The only tools you would need for this are the BLAST executables (found here) and the databases (info on how to download can be found here).
I should note that you don't really need any Python to get this to work, only some shell scripting. Although you can run these commands through a Python wrapper if you really want to.
Good luck!

So that's like 60M BLAST queries? How long would that take?

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


text:
Extracting features or gene from PCA after calculating PCA for downstream analysis
The quote below is from this paper:

We performed principal component analysis (PCA) of low-coverage sequencing data to identify genes explaining variation across cells. PCA separated the cells into groups corresponding to the source populations (Fig. 2c and Supplementary Figs. 3–5), and the genes distinguishing each group reflected the biological properties of the cell types (Supplementary Fig. 5 and Supplementary Table 3). PCA of low- and high-coverage sequencing data revealed a remarkably similar graphical distribution of analyzed cells, and the majority (78%) of the top 500 genes determined by PCA were shared between PCA performed on low- and high-coverage data (Supplementary Figs. 4 and 6 and Supplementary Table 4).

How are they performing PCA and then finding out the components which contribute to the strongest PC and then taking out the genes for downstream analysis? Let's say there are genes involved in PC1 and PC2 which define a certain cell type and that distinguish it. I think for a biologist, like me, I just want to get the genes or gene list that are getting involved in determining lineage or cell type. How can I do this?
As a to test what they might be doing i did this
test<- read.csv('NON_CODING.csv',header = T,row.names = 1)

# preform PCA
pca = prcomp(t(test), center=TRUE, scale=TRUE)

then after calculating PCA i do see this in my console when I type pca
> pca
Standard deviations (1, .., p=16):
 [1] 3.826851e+01 2.080405e+01 1.568739e+01 1.349256e+01 1.119348e+01 9.980547e+00 8.365034e+00
 [8] 8.098841e+00 7.519507e+00 6.184505e+00 5.880260e+00 5.139609e+00 4.851091e+00 4.335606e+00
[15] 3.870918e+00 3.173450e-14

Rotation (n x k) = (2899 x 16):
                                       PC1           PC2           PC3           PC4
5S_rRNA                      -1.090574e-02 -2.412665e-03  1.637689e-02 -3.603865e-02
AB019441.29                  -1.928250e-02  1.821083e-03 -9.713724e-03 -1.978737e-03
ABBA01017803.1               -1.823266e-02 -3.727144e-03 -9.790131e-04 -3.937062e-02
ABC14-1080714F14.1            2.438024e-02  3.816019e-03  1.657784e-04 -7.141480e-03
ABC7-481722F1.1               2.467403e-02  2.873432e-03  3.781153e-03  3.790824e-03
AC000036.4                   -1.066777e-02  2.838305e-02  2.642673e-02 -7.998608e-03
AC000089.3                   -1.249026e-02 -7.870864e-04 -2.569602e-02 -3.165073e-03
AC000120.7                   -1.313696e-02  5.923245e-03  1.633255e-02  4.077602e-02
AC000123.4                    2.415996e-02  1.029732e-02  8.930792e-03 -1.118660e-02
AC000403.4                    8.692835e-03 -3.060170e-02  3.746839e-02 -6.968971e-03
AC002064.4                   -1.799590e-02  2.147160e-02 -2.303742e-02  7.647937e-03

Now how do i find which are the list or set of genes that gives me PC1 or PC2 so on
Any suggestion or help how to get the genes that make most difference between two major component that would be very helpful

The help page explains that there is a component which:

rotation
  the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.

So using using the rotation one can use it to calculate the importance of each feature in the input:
C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
pZ <- prcomp(Z)
genesImportance <- Z %*% pZ$rotation
genesImportance2 <- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"

This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data.

Im using prcomp function would that work as you mentioned princomp ?

Z%*%pZ$rotation should be the same as pZ$x (that’s in fact how pZ$x is calculated by prcomp).

Totally true, I didn't went further in the help page, but with rotation one can predict the position of new samples not included in the original data. I'll edit it.

@KonradRudolph I'm still not getting where is this coming from  italic pZ$x   im little ignorant as I see examples and do it so i have to do a back referencing  to in R methods and tutorials , is it in my data frame as well ?with context to my data frame are you suggesting this  italic"pca$x"  and what is that  Z  actually

@krushnachChandra Z is the initial data for the PCA. pZ is the output of prcomp. See the edit, I might not explain it clearly enough

okay i will try to find it out

@Llopis i tried that i get error like this Error in Z %*% pZ$rotation :
requires numeric/complex matrix/vector arguments

@krushnachChandra See the new edit. Could you explain what code are you using and what is Z?


text:
Exception : [Failed to read 8 bytes from input stream! Read 0] salmon quant was invoked improperly
I am trying to run salmon and it keeps giving me 2 java exceptions:
java.io.EOFException
    at java.util.zip.GZIPInputStream.readUByte(GZIPInputStream.java:268)
    at java.util.zip.GZIPInputStream.readUShort(GZIPInputStream.java:258)
    at java.util.zip.GZIPInputStream.readHeader(GZIPInputStream.java:164)
    at java.util.zip.GZIPInputStream.<init>(GZIPInputStream.java:79)
    at java.util.zip.GZIPInputStream.<init>(GZIPInputStream.java:91)
    at uk.ac.babraham.FastQC.Utilities.MultiMemberGZIPInputStream.<init>(MultiMemberGZIPInputStream.java:37)

and
java.lang.ArrayIndexOutOfBoundsException: -1
    at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.calculateDistribution(SequenceLengthDistribution.java:100)
    at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.raisesError(SequenceLengthDistribution.java:184)
    at uk.ac.babraham.FastQC.Report.HTMLReportArchive.startDocument(HTMLReportArchive.java:336)
    at uk.ac.babraham.FastQC.Report.HTMLReportArchive.<init>(HTMLReportArchive.java:84)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.analysisComplete(OfflineRunner.java:155)
    at uk.ac.babraham.FastQC.Analysis.AnalysisRunner.run(AnalysisRunner.java:110)
    at java.lang.Thread.run(Thread.java:745)

    at uk.ac.babraham.FastQC.Sequence.FastQFile.<init>(FastQFile.java:80)
    at uk.ac.babraham.FastQC.Sequence.SequenceFactory.getSequenceFile(SequenceFactory.java:106)
    at uk.ac.babraham.FastQC.Sequence.SequenceFactory.getSequenceFile(SequenceFactory.java:62)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.processFile(OfflineRunner.java:129)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.<init>(OfflineRunner.java:102)
    at uk.ac.babraham.FastQC.FastQCApplication.main(FastQCApplication.java:316)

I installed salmon using Anacondas' conda install -c bioconda salmon and all other necessary packages in the same way. I thought that it could be a problem with fastqc, so I uninstalled it and then installed it manually (through the fastqc.zip file), but the output remained the same.
I am using Java8:
java -version
openjdk version "1.8.0_92"
OpenJDK Runtime Environment (Zulu 8.15.0.1-linux64) (build 1.8.0_92-b15)
OpenJDK 64-Bit Server VM (Zulu 8.15.0.1-linux64) (build 25.92-b15, mixed mode)

I am running the script on a cluster with SLURM. Locally on my laptop the script is running without any issues. I am really stuck and do not know what to try, so any suggestions would be greatly appreciated.
Update
Another issue that I just found is the following (the last command that salmon is trying to execute):
Exception : [Failed to read 8 bytes from input stream! Read 0]
salmon quant was invoked improperly.

Then I searched online for the issue and found the following thread:
https://github.com/COMBINE-lab/salmon/issues/129
However, updating salmon to the newest 0.9.1 version did not solve the issue.
Edit
The command I am running is the following one:
python "$path/salmon_rna_seq.py" -L 'IU' -E "$path/jyvo_experiment-
metadata.yaml" -T 'gencode_mouse_m13' -C "$path/salmon_rna_seq.yaml" -I
"$path_to_fastq/$name-lec_1_r1.fq.gz" "$path_to_fastq/$name-
lec_1_r2.fq.gz" -O "$path/res/$name-lec_1" -S "$name-lec_1" -P '6' -X 'True'

The command is run by Pypiper: https://github.com/epigen/pypiper
Inside of salmon_rna_seq.py I am also running commands to treat the data, the salmon part looks the following way:
def salmon(self):
    cmd = "salmon quant -i " + self.indexed_transcriptome
    cmd += " -l " + self.lib_type

    if self.trim:
        cmd += " -1 " + self.sample_files['r1']['fastq']['trimmed_paired']
        cmd += " -2 " + self.sample_files['r2']['fastq']['trimmed_paired']
    else:
        cmd += " -1 " + self.sample_files['r1']['fastq']['chastity']
        cmd += " -2 " + self.sample_files['r2']['fastq']['chastity']

    cmd += " -o " + self.sample_dirs['salmon']['base']
    cmd += " --numBootstraps=30"
self.pipe_manager.run(cmd, self.sample_files['salmon'])

Commands that the Pypiper runs are the following ones (_commands.sh Pypiper file):
zcat -c $PATH/res/old-lec_1/fastq/merged/old-lec_1_r1.fq.gz | grep -A 3
'^@.*[^:]*:N:[^:]*:' | grep -v '^\-\-$' | gzip -c > $PATH/res/old-lec_1/fastq
/chastity/old-lec_1_r1.fq.gz

zcat -c $PATH/res/old-lec_1/fastq/merged/old-lec_1_r2.fq.gz | grep -A 3
'^@.*[^:]*:N:[^:]*:' | grep -v '^\-\-$' | gzip -c > $PATH/res/old-lec_1/fastq
/chastity/old-lec_1_r2.fq.gz

java -jar $HOME/anaconda2/share/trimmomatic-0.36-5
/trimmomatic.jar PE -threads 6 -phred33 $PATH/res/old-lec_1/fastq/chastity/old-
lec_1_r1.fq.gz $PATH/res/old-lec_1/fastq/chastity/old-lec_1_r2.fq.gz $PATH/res/old-
lec_1/fastq/trimmed/old-lec_1_r1.fq.gz $PATH/res/old-lec_1/fastq/trimmed/unpaired
/old-lec_1_r1.fq.gz $PATH/res/old-lec_1/fastq/trimmed/old-lec_1_r2.fq.gz $PATH/res/old-
lec_1/fastq/trimmed/unpaired/old-lec_1_r2.fq.gz ILLUMINACLIP:$HOME/anaconda2/share/trimmomatic-0.36-5/adapters
/Small_List_Of_Adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15
MINLEN:36

fastqc $PATH/res/old-lec_1/fastq/merged/old-lec_1_r1.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/merged

fastqc $PATH/res/old-lec_1/fastq/chastity/old-lec_1_r1.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/chastity

fastqc $PATH/res/old-lec_1/fastq/merged/old-lec_1_r2.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/merged

fastqc $PATH/res/old-lec_1/fastq/chastity/old-lec_1_r2.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/chastity

fastqc $PATH/res/old-lec_1/fastq/trimmed/old-lec_1_r1.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/trimmed

fastqc $PATH/res/old-lec_1/fastq/trimmed/old-lec_1_r2.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/trimmed

salmon quant -i $HOME/Genomes/mouse/salmon
/gencode.vM13.transcripts.fa.sidx -l IU -1 $PATH/res/old-lec_1/fastq/trimmed
/old-lec_1_r1.fq.gz -2 $PATH/res/old-lec_1/fastq/trimmed/old-
lec_1_r2.fq.gz -o $PATH/res/old-lec_1/quantification --numBootstraps=30

Edit
After running manually gzip -t on a file:
-bash-4.2$ gzip -t old-lec_1_r1.fq.gz
gzip: old-lec_1_r1.fq.gz: unexpected end of file

So, I guess the file is corrupted and the issue is with some library that is generating the file. However, reinstalling manually other libraries do not help.

Can you run gzip -t on the file(s) to see if it/they are corrupt? That'd be the interpretation of the error messages you're seeing. Make sure to run this on the cluster, in case there's some sort of "modified in transit" networking error going on.

This would be my guess.  Also, however, I'm not sure how Java is involved (since Salmon is written completely in C++).

It looks like FastQC is being invoked somewhere. Providing the precise commands invoked will help us troubleshoot the problem.

Yes, please [edit] your question and show the exact commands you are using. It's hard to debug something if we don't know how to replicate it.

Salmon is a red herring. The error message is clearly unrelated to Salmon.

Oh god. I swear I didn’t notice the pun in my previous comment.

+1 Konrad, for the unintentional pun ;P.

Where does old-lec_1_r1.fq.gz come from?  That is the issue.

FastQC is involved in generating it mainly

FastQC has absolutely nothing to do with generating the file, it just produces some QC metrics. It's unclear (A) why you're running grep on the files and (B) which of the multiple copies of old-lec_1_r1.fq.gz is the first that's corrupt.

The issue was caused by the corrupted input .fastq files that were damaged somehow upon uploading them from the local machine. We figured it out by md5sum command output comparison. We have not yet managed to actually run all of the script still, it is failing, but for another reason now.


text:
t-sne code and labelling issue for rna seq data
I'm trying to run t-sne on my coding as well as non coding genes. Even though I'm getting output I'm not confident about my code (The output it's difficult to interpret).
I have a data set of total 16 cell type and about 3000 genes.
This is how my data is first column is gene rest is sample:
gene    H1  H2  H3  H4  M2  M3  M4  C1  C2  C3  C4  G2  G3  G4  Gran1   Gran2

The code used:
library(Rtsne)
file1<- read.csv('UNIQUE_WGCNA_FILE.csv',header = T,row.names = 1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=50, theta=0.5, dims=3)

cols <- rainbow(10)
plot(tsne$Y, t='n')
text(tsne$Y, labels=as.numeric(trn[,1]), col=cols[trn[,1]])

require(rgl)
plot3d(tsne$Y, col=cols[trn[,1]])
legend3d("topright", legend = '0':'5', pch = 16, col = rainbow(10))
rgl.postscript("Merge.pdf", "pdf")

I'm getting output but I don't know how to label so its becoming difficult.
Any suggestion or help how to label would be appreciated.

What kind of labels are you looking for?

I want to label them according to my sample such as I have H1  H2  H3  H4  M2  M3  M4  C1  C2  C3  C4  G2  G3  G4  Gran1   Gran2

@krushnachChandra I've edited the question, I hope I have kept the meaning. Could you add the image of the T-SNE without the labels? Also I would like to know how have you tried to solve this (I don't want to suggest things you already tried).

@Llopis first of all my code is correct? or not let me know as you know the dimension of my input , so the code is correct for my input

@Llopis i have added my image and i saw the Rtsne manual i can't find ways to get the column name and label it in the plot ..so i couldn't do it..any help would be appreciated

Although a bit late, here is the solution for this:
1) Create color/column levels:
levels = unlist(lapply(colnames(trn[,1:16]), function(x)
    {

        if(x == "H1") return("H1")
        else if(x == "H2") return("H2")

}))

2)
plot(tsne$Y, col=factor(levels))

This should color columns H1 and H2 differently
Edit:
This is a sample complete code:
library(Rtsne)
file1<- read.csv('filename.csv',header = T, row.names=1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=5, theta=0.5, dims=3)

collevels = unlist(lapply(colnames(trn), function(x)
{
  if(grepl("H1",x))
    return("blue")
  else if(grepl("H2",x))
    return("red")

  # etc ... for all the columns you want to color
}
))

colorColumns = factor(collevels)
require(rgl)
plot3d(tsne$Y[,1], tsne$Y[,2], tsne$Y[,3], col=colorColumns)

thank you for the code i will give it a try, but do you have full tnse code if yes i would be glad if you can post in the answer i will give a try

It seems that tsne$Y is the matrix and the names are just taken from there. So, you could replace the names of the columns of that matrix and then it would display the names instead of the numbers.
But, you are using the function plot3d from the rgl package, which has as arguments xlab, ylab, and zlab, so you just need to give a string to those arguments to name the axes.

i want to label my sample after the plot the axis isn't a problem ,because it's difficult to identify which point coming from which sample ..

Do you think that adding the labels of 3000 genes will help to visualize the plot? How will this help you?

nope i just want to lalel the sample name....not genes....I mean based as I have  samples as these H1  H2  H3  H4  M2  M3  M4  C1  C2  C3  C4  G2  G3  G4  Gran1   Gran2 ,the H1-H4 hsc, then M2-M4 monocyte , C1-C4 CMP , G2-G4 GMP and then Gran1 and Gran2 are Granulocyte ..I want to label them in terms of my sample name ..Particular color corresponding to particular group of sample  ,I hope now you can understand my issue

Do you realize that you have more than 30 points in your t-SNE? No, I can't understand your issue, please edit the question and explain further what do you want

okay so I will rephrase it again so just like in PCA you can label a group of sample that cluster together lets say same lineage , I would like to do the same here too.

So, could you edit your post including an image similar to the desired output? I don't grasp what do mean.

okay will do that


text:
Finding a single open reading frame with ribosomal binding site, using Biopython
I'm given a Fasta file, containing a large DNA(over 115,000 long) sequence, and I am tasked with finding a single large open reading frame contained within the DNA sequence using Biopython.
I'm aware this has been asked before consistently however i can't find anything for actually finding the source of my ORF, but apologies for any overlap.
Now from other sources and the Biopython cookbook I've translated my sequence and found six open reading frames (three for each strand) and their positions within the sequence;
def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    seq_len = len(seq)
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            trans = str(nuc[frame:].translate(trans_table))
            trans_len = len(trans)
            aa_start = 0
            aa_end = 0
            while aa_start < trans_len:
                aa_end = trans.find("*", aa_start)
                if aa_end == -1:
                    aa_end = trans_len
                if aa_end-aa_start >= min_protein_length:
                    if strand == 1:
                        start = frame+aa_start*3
                        end = min(seq_len,frame+aa_end*3+3)
                    else:
                        start = seq_len-frame-aa_end*3-3
                        end = seq_len-frame-aa_start*3
                    answer.append((start, end, strand,
                                   trans[aa_start:aa_end]))
                aa_start = aa_end+1
    answer.sort()
    return answer

orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)

for start, end, strand, pro in orf_list:
print("%s...%s - length %i, strand %i, %i:%i" \
      % (pro[:30], pro[-3:], len(pro), strand, start, end))

OUTPUT
 TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427

However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of AAGGAGGTG and occurring between 6 to 9 nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?
My specification says to find the single open reading frame however since none of these begin the the M start codon i'm a bit lost how they can be so, I have no information on where my sequence is from (something I have to find out) so using the standard translation table 1.

You should start by finding all instances of AAGGAGGTG and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.
There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading frames. You seem to be starting the reading frame at any non-stop codon instead of starting from AUGs. This may well make biological sense in some contexts (there are proteins that don't start with M) but is probably beyond the scope of what sounds like a simple assignment.

Thanks for the heads up, yeah I think I can safely assume that yeah the start codon will be an M. Apologies for the follow up question, is there anywhere you can point me in the right direction to in terms of achieving this? I'm struggling to find anything related.

@maltedg sorry, I just saw your comment, but I don't know what you mean. This is pointing you in the right direction. As I said, start by finding AAGGAGGTG, and then look for ORFs that start after it. If you can assume that they'll always start with an AUG, then look for AUGs close to the AAGGAGGTG and extend the reading frame from there.

From what I can tell, the program output does not indicate whether the open reading frames (ORFs) are preceded by ribosome binding sites (RBSs). I'm not sure what parameters are used when establishing ORFs. To verify your RBS requirement, index the input DNA (using given ORF locations from biopython output, pay attention to strands, which are designated by 1/-1) and search upstream for a sequence similar to your consensus. You'd have to pre-define your criteria for similarity to consensus in order to find it programmatically (or you can do it visually to better avoid false negatives). Perhaps try a search that permits ambiguous/non-matching third nucleotides in codons (like a regular expression such as "AA.GA.GT." where "." can be any nucleotide). Biologically, third nucleotide changes are least probable to drastically alter the translated residue (more likely to encode amino acid with similar chemical properties).


text:
Which tools are available to call somatic mutations on non-coding regions from whole genome sequencing data?
I have a set of whole genome sequencing samples, some of which have matched normal while some not. I want to call somatic mutations on non-coding regions.
I was looking for GATK best practices (like this) but it seems it is not yet ready for production. Which alternative tools do you suggest?

Have you looked at this question? https://bioinformatics.stackexchange.com/questions/2744/somatic-tumor-only-variant-calling

@Bioathlete thanks, I actually didn't find that question before. However, I found that question too specific on tumor-only samples. I also have normals so maybe there are better tools?

Those will have to be treated as separate analysis.  You can use Mutect or VarScan for the matched samples and the tools listed in the other question for the unmatched samples however that means that you won't be able to compare the results easily.  If the inter-sample comparison is important you will have to run all samples as unmatched.  Personally I have used LoFreq for somatic variant detection with some success.  It is very sensitive so the false positive rate can be high depending on how noisy your data is.


text:
Get gene annotation type from gene name
I am having RNA-seq data set for C.Elegans. I do not know which annotation scientists used for naming the genes of it, but I need to search somehow for the ^ERCC and ^mt- genes like specified in the tutorial:
https://www.bioconductor.org/help/workflows/simpleSingleCell/
If I just implement grepl search for the character array of gene names that I have, I am always getting all values FALSE as if there are no spike-ins and mitochindrial genes.
Any suggestion would be greatly appreciated.
Sample gene name from the dataset: WBGene00023193

The tutorial you linked to does not mention ERCC or mt genes.
Before running grepl or any other search, you should first check your inputs. It looks like you have WormBase gene IDs if all your genes look like WBGene00023193. Therefore, your question is actually how to convert WormBase gene IDs to gene symbols.
Here is one potential solution: https://github.com/WormBase/website/issues/4297
After you convert, you have to make sure that you have ERCC spike-ins in your gene table (probably not) and you mitochondrial genes start with "mt-".

Sorry about that, the link was wrong. I updated it


text:
Downloading URL using AWK for fixed fields
I have a text file as follows
a   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz
b   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz
c   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz

The file is tab delimited and has 8 columns with 8th column being the url link  I need to download the files one by one using wget and changing the name of the downloaded file with -O option in wget. The argument I want to give to the -O option is combination col1 and col2 separated by _ so that for first line I get the downloaded file as a_k.fastq.gz and so on. For example, the first line of the file above should give:
wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz

Following the post from here I realized looping through the file is slow and sometimes difficult. How can I download these files using wget with awk or cut as the tool for fixed fields

Can you post up an example wget command?

wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz

how can I put it into code in the comment section??

edit the original post

@Llopis we've decided that basic things that bioinformaticians need to do in their daily work are on topic even when the issue could be expressed as a general text parsing question.

Ok, I misjudged and thought it was about downloading data.

awk '{print "wget " $9 " -O " $1 "_" $2 ".fastq.gz"}' tmp

will give you
wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz -O b_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz -O c_k.fastq.gz

If you add a | bash to the command it will then execute all the wget commands.

This works great. Can you tell what |bash is doing ??

@user3138373 the output of the awk command is a set of shell commands, essentially a mini script. By piping it to bash with | bash you're telling it to execute each input line as a command. So it will run the wget calls and download the sequences.

A variation of Bioathlete’s answer that doesn't suffer from escaping issues:
awk 'BEGIN{ FS="\t" } { printf("%s\t%s_%s.fastq.gz\n", $8, $1, $2) }' input |
while IFS=$'\t' read -r url filename; do
    wget -O "$filename" -- "$url"
done

Or entirely with shell parsing:
while IFS=$'\t' read -r c1 c2 c3 c4 c5 c6 c7 url; do
    wget -O "${c1}_${c2}.fastq.gz" -- "$url"
done < input

Or entirely without shell parsing but with GNU xargs:
awk 'BEGIN{ FS="\t" } { printf("%s_%s.fastq.gz\n%s\n", $1, $2, $8) }' input |
xargs -r -d '\n' -n 2 -- wget -O

… you could even run multiple wget instances in parallel from GNU xargs with the -P option, though you should set the -q option on wget to disable progress bars overwriting each other on the same terminal. Then again it would be better to use a real command-line download manager like Aria or Axel for that task.


text:
Where can I get the population allele frequency vcf file?
I want to use GATK to estimate cross-sample contamination for Whole Genome Sequencing data.
The specific tool is ContEst and it is run with:
java
 -jar GenomeAnalysisTK.jar \
 -T ContEst \
 -R reference.fasta \
 -I:eval tumor.bam \
 -I:genotype normal.bam \
 --popfile populationAlleleFrequencies.vcf \
 -L populationSites.interval_list
 [-L targets.interval_list] \
 -isr INTERSECTION \
 -o output.txt

The --popfile option requires a vcf file representing population allele frequencies. Where can I get it? Is it the one I can get from dbSNP?
I'm using the human_g1k_v37.fasta reference, a version of hg19, "for the phase1 analysis we mapped to GRCh37. Our fasta file which can be found here called human_g1k_v37.fasta.gz, it contains the autosomes, X, Y and MT but no haplotype sequence or EBV".

On the GATK forum they've recommended the population stratified VCF file for this purpose.


text:
Determine if a gene is mitochondrial or not
I need to determine whether a gene is mitochondrial or not for C. Elegans automatically from its name by running a regular expression on a dozen of thousands of gene names. Currently I am having gene symbols like:
clec-190, clec-189, 21ur-8912, nduo-1, 21ur-6539, 21ur-13096, nduo-6, 21ur-14241, ndfl-4

I tried running grepl("^mt-", symbols) and it found nothing. I am sure that there are some mitochondrial genes just from their descriptions. I would assume that first I need to convert gene symbols to some other symbols to do that but I am not sure here.
Basically, I am trying to do it in a similar way to the tutorial of single cell RNAseq.

Please [edit] your question and clarify what species and what sort of gene identifier you will be using. In any case, the chances of this being the type of problem you can solve with regular expressions are vanishingly small. Is this something you will need to do once or many times? It would be much simpler to query online databases for information on the gene rather than try to parse the names.

Ok, I added the species and a bit of clarification

I also have Worm Base identifiers for the genes if that can help

I know that the gtf files provided by illumina for C elegans use Wormbase identifiers ("WBGene*"): http://jp.support.illumina.com/sequencing/sequencing_software/igenome.html. The gtf files contain the chromosome, so with tools like grep and awk, you can easily know wheter a given identifier is from the mitochondrial chromosome.

You can't determine if a gene is mitochondrial purely from its name. The simplest route would be to use the wormbase biomart and either download the gene names and their associated chromosomes or use the biomaRt package from bioconductor to do the same thing. You can upload a list of the genes you're interested in querying in either case (or just download all of them, and filter).
From your list, only nduo-1, nduo-6, and ndfl-4 are mitochondrial.


text:
Intersecting two different files with one "master" file based on different columns
I have the following sets of data:
file1:
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    16855676    16855696    CTGTCCAGCAGAGGGCGGTG

file2:
1   15776231    2   5008    G:5002  A:6
1   15776239    3   5008    C:3358  A:14    G:1636

file3:
1   16855677    2   5008    A:5003  C:5
1   16855689    3   5008    G:4957  A:41    T:10

Most lines have 6 columns but a few lines also have 7 columns which have data that I want to preserve.
Can a script be written which does the following:
if $2 of file2 comes in between $2 and $3 of file1 (between the interval of $2 and $3) and $2 of file3 comes in between $5 and $6 of file1 (between the interval of $5 and $6), then the output would be $1,$2,$3,$4 from file1 and $2,$3,$4,$5,$6,$7 from file2 and $5,$6,$7 from file1 and $2,$3,$4,$5,$6,$7 from file 3, all in one line. So something like this:
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776239    3   5008    C:3358  A:14    G:1636  16855676    16855696    CTGTCCAGCAGAGGGCGGTG    1   16855689    3   5008    G:4957  A:41    T:10

Something like bedtools does but I am not able to use bedtools one this since I don't have two columns in file2 and file3. The intersection is based on $2  in both the files, and all the other columns after it are just copied
(There could be many lines from file2 and file3 intersecting in the same line in file1,so that line could be repeated based on the number of intersections)

Can't you just munge file2 to add an end coordinate and then use bedtools intersect? For example awk '{$3 = $2+1 FS $3}' file2 > file2.modified, assuming the coordinates are 0-based.

Can you please clarify how this is related to your previous question?

@heathobrien it's an extension. When the OP asked the previous question, he hadn't realized there were certain issues in the input files. So I suggested he ask a new question showing the files he currently has after fixing these issue. This is now about getting information from three separate files and two formats (number of fields) per file. So quite a different issue.

An awk solution
awk '
    FILENAME == ARGV[1] {f2[$2] = $2 OFS $3 OFS $4 OFS $5 OFS $6 OFS $7; next}
    FILENAME == ARGV[2] {f3[$2] = $2 OFS $3 OFS $4 OFS $5 OFS $6 OFS $7; next}
    {
        for (key2 in f2) {
            for (key3 in f3) {
                if (($2 <= key2 && key2 <= $3) && ($5 <= key3 && key3 <= $6)) {
                    print $1,$2,$3,$4, f2[key2], $5,$6,$7, f3[key3]
                }
            }
        }
    }
' file2 file3 file1

outputs
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776231 2 5008 G:5002 A:6  16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855677 2 5008 A:5003 C:5
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776231 2 5008 G:5002 A:6  16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855689 3 5008 G:4957 A:41 T:10
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776239 3 5008 C:3358 A:14 G:1636 16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855677 2 5008 A:5003 C:5
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776239 3 5008 C:3358 A:14 G:1636 16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855689 3 5008 G:4957 A:41 T:10

For your sample data, all combinations of records in file2 and file3 meet the criteria for the file1 record.


text:
How to filter ribosomal RNA from scRNA-seq data
I want to filter out ribosomal RNA from scRNA-seq data (downloaded from here). Is there a list of known ribosomal RNA?
The only solution I found is SortMeRNA, however it works with raw sequencing data afaik, while I already have a matrix with transcript counts for each cell. I searched for a comprehensive list of rRNAs but I didn't find any.

What GTF file did you use to generate the counts? For some species only 5S (and maybe MT) rRNA is left, for others you might have one or more copies of the various 48S components.

I actually don't have GTF files. I downloaded the data provided by the Tabula Muris Consortium (https://figshare.com/projects/Tabula_Muris_Transcriptomic_characterization_of_20_organs_and_tissues_from_Mus_musculus_at_single_cell_resolution/27733)

The rRNA genes in that dataset are Rn45s and Rn4.5s.
BTW, you have gene counts, not transcript counts.

Thanks. What do you mean by gene counts instead of transcript counts? That the counts are already collapsed by gene?

Also, how did you find the list of rRNA genes in such dataset?

@gc5 The counts were likely only ever done by genes to begin with. rRNA genes typically have names starting with "Rn" in mouse.

Ok. Should I need to find ribosomal genes in human do you know any database which I can use?

I think they have similar names in humans.

In the paper mentioned, we used the ScaleData function in Seurat to regress out the number of reads, Rn45s abundance, and percent ribosomal gene transcripts. Ribosomal genes were found with the regular expression ^Rp[sl][[:digit:]].
tiss <- ScaleData(object = tiss, vars.to.regress = c("nReads", "percent.ribo","Rn45s"))

Here's a fuller notebook, and we'll have a better organized repository soon.

That's great that the workflow is available. I am just curious why did you decide to use ribosomal genes instead of the mitochondrial genes, which is suggested by Seurat?

Hello! Thanks for asking. In the final version of the paper (this question cites the preprint), we ultimately decided not to normalize on either ribosomal or mitochondrial genes. We found that this normalization added differences where there shouldn't be any, e.g. between cells of the same type in the Bladder.

Thanks for clarifying! I was wondering why I could not find any mention of ribosomal genes anywhere in the methods (in the final paper).


text:
Difference between de novo transcriptome assembly methods
I have been looking around (including read the original papers) to understand what is essentially the difference between StringTie in non-reference based mode (de novo) and Trinity de novo assembly. I understand that in the genome-guided nature of StringTie, we are supported by the genome information (i.e. taking into account the mapped reads). In addition to that, we also have an option to either run using the gene annotation file or not.
Simply put, am I correct if I say that Trinity is purely looking at the overlapping reads, then assemble them, whereas StringTie (non-reference based) is looking at the mapped reads that are either close or overlapping to each other? Is there any other thing I'm missing here (I understand it's not that simple, but I'm trying to make it as simple/intuitive as possible)?
Also if that's the case, would it be safe to say that generally non-reference based StringTie is preferred over Trinity and reference based StringTie (as gene annotation can lead to some "bias", hence preventing discovery of novel TSS, exon length or even make unnecessary duplicates with small differences)?

Your intuition is correct. stringTie is just looking at clumps of alignments and how they might relate to each other (either due to spliced alignments or proximity). Trinity is doing the more computationally difficult task of finding parts of reads that overlap with each other and trying to link them together into longer sequences.
Whether it makes sense to use something like stringTie or not depends mostly on the quality of the reference sequence. If you're working on an organism with a high-quality reference sequence (pretty much anything with a lot of papers using it) then there's no reason to bother with Trinity. If your reference genome additionally has a high-quality annotation (human, mouse, etc.) then you should use it, as anything "novel" that the program is likely to output otherwise is probably just an artefact. If you either don't have an annotation or it's known to be highly problematic then it makes sense to omit it when using stringTie.
Without seeing evidence to the contrary, I don't buy the argument that using a good annotation is actually biasing the output in a problematic way. For well-studied organisms, anything reported as "novel" these days likely doesn't exist, or is otherwise so fleeting that it's biologically irrelevant. The arguments of bias sound like something that was come up with once someone didn't get as many "novel" results as they were hoping.

Thanks for your answer. About the point we will less likely find novel things in organism with high-quality annotation, I was wondering if that's always the case as there are many recent publications highlighting novel elements, even in human. I guess you are saying that generally it's sufficient, but in the case of some specific condition or localisation, finding novel elements is still relevant?

Finding a few of them, sure, but unless the system is being severely disturbed one would only expect a handful of real novel findings (those reporting more are likely primarily false-positives at this point).

It's really a misnomer to call StringTie's non-reference based mode 'de-novo.' It's still using the reference genome sequence to guide the transcript assembly, it's just not using the reference annotation. Trinity is truly de-novo in that it it assembles the transcript from the overlap of the reads without mapping them to a reference genome sequence.
If you look at the precision numbers from the StringTie paper, you'll see that it maxes out around 80% under ideal circumstances with simulated data. This means that if you're working with a well annotated genome where most transcripts have already been discovered, the vast majority of the 'novel' transcripts you find will likely turn out to be artifacts.

Thanks for the answer. I agree for the first point. For the second point, I partially agree as I think we can always find some novel transcripts in the case of uncovered condition/localisation (as my comment to Devon's answer). I hope "vast majority" here means that there could still be a few out of thousands that are indeed biologically relevant. ;)

Happy to help. It's certainly possible that some novel transcripts will be real, especially if your experimental setup is novel, but validation is likely to be a major bottleneck. In my opinion, short reads are really only suitable for quantification of known transcripts. Transcript discovery is much more effective with long read sequencing technologies like PACBIO Iso-Seq or Oxford Nanopre.


text:
What are RIKEN genes?
I am using data from the Tabula Muris Consortium.
These are gene counts from scRNA-seq of mouse cells.
There are some specific genes with name ending with suffix 'Rik' (e.g. 0610005C13Rik or 0610007C21Rik). Searching online I found that they are from RIKEN project.
However, they are different in some respects (e.g. antisense lncRNA gene, or mouse 0610007C21RIK Adenovirus).
I want to exclude them from downstream analysis, but it is unclear to me the main purpose of this collection. For which kind of analysis are they used?

My understanding is that these are cDNA sequences from the RIKEN project that don't yet have another gene name. The project carried out a comprehensive sequencing of the mouse genome, and as a result found a lot of cDNA transcripts that were previously unknown. Now they have some annotation (which can be varying degrees of functional annotation), and an accepted official identifier to allow people to talk about these transcripts when discovered in other studies.
It'd be interesting to see if any novel transcripts have been found from the Tabula Muris project. They would also need their own identifier, possibly leading to similar intrigue about their history.

cDNA is generated from RNA so they most probably sequenced RNA (the transcriptome) instead of the genome itself.

To extend a bit more to @gringer answer. RIKEN is a large Japanese research institute, that participates in different scientific collaborations.
One of these is the FANTOM collaboration (Functional ANnotation Of the Mammalian genome), which lately have resulted in the release of FANTOM CAT, a browsing tool for exploring lncRNA relative to associated Genes, Traits and Ontologies.
Many of the transcripts and genes discovered during this work, are not yet named, and are therefore called CAT-genes, following numbering comparable to ensembl (eg. CATG00000118335.1).
I'm unsure how RIK genes have been discovered, as I have not personally encountered these genes yet, but in case you likewise encounter more genes discovered in collaboration with RIKEN, you now also know of CAT genes.


text:
Normalization using parallel R script
I need to implement a parallel version of an R script because on one core it takes forever to execute. The function that takes the most of time is scran::quickCluster() applied to an integer sparse matrix (UMI_count) that is 26000x66000. The code is presented below:
library(scran)
library(scater)
sce <- newSCESet(countData=UMI_count)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, cluster=clusters)

I did not find much on parallel versions of quickCluster() function, so will it actually execute in parallel if I set the number of cores to more than one? I am trying to do normalization of heterogeneous dataset with multiple cell types, that is why I am using quickCluster().
And I just need to speed it up to be able to compute it.

Do you really have to use a sparse matrix? A matrix of that size with 64-bit integers is still <20GB and not have any of the annoyances of a sparse matrix.

But will it speed up quickCluster() function in any way? I can convert it to regular matrix on the cluster (that gives up to 120Gb of RAM).

I suspect it will, it's generally easier to work with normal matrices.

The scran::quickCluster method has a BPPARAM (i.e., BiocParallelParam), so if one provides something like
clusters <- quickCluster(sce, BPPARAM=BiocParallel::MulticoreParam())

it will use all available cores (a specific number of cores can also be given). See the intro doc to BiocParallel for info on other options.


text:
Order of batch effects removal, data imputation and library size normalization in scRNA-seq data
I am preprocessing scRNA-seq data. What is the best practice in use to run both ComBat for batch effects removal, data imputation (to mitigate dropout) and library size normalization?
I thought that library size should be run first, since it is per-cell normalization, then ComBat batch effects removal. On the original paper - Johnson et al. (2007) - it is stated that:

We assume that the data have been normalized and expression values have been estimated for all genes and samples.

However, I want to apply it to scRNA-seq data. Does this statement still hold? Additionally, I plan to apply imputation (e.g. with MAGIC) in the end. Is there any problem you can spot?
Update
I attach the PCA regarding an example Mus Musculus dataset in which different colors represent different mice. It seems clear to me that the first two principal components are affected by batches (mouse id).

Update 2
I rerun the PCA on raw counts data (the first PCA was on log-transformed data) and I obtain a different description of the dataset, in which batch effects seem not to be prevalent.

From what I can tell, MAGIC should be run on raw data, so that would be the first step.

@burger MAGIC normalizes data before imputation, so it should be run at least after library size normalization. My concern is that using MAGIC before ComBat will amplify batch effects. Reading the paper I could not find any reference to batch effect removal.

The advice I got was that the best would be to adjust for batch effect instead of removing them. Did you tried adjusting for your batch effects? How big is your batch effect ? (Do the PCA or MDS or dendograms show a clear distinction by your batch effect (or several batchs) ?)

@Llopis yes, actually for batch effects removal I meant adjusting for batch effect with ComBat, is it what you meant?

No, comBat doesn't adjust for batch effect it "removes" it (Despite the first line of the help page). From the (same) help page: "Users are returned an expression matrix that has been corrected for batch effects"; it modifies the data to "adjust" for it, instead of adding/calculating a factor to be taken into account for later steps. The later can be done in limma, DESeq2 and other packages but it is not the same adjusting than removing.

@Llopis ok thanks, I didn't know this distinction. However, my PCA shows a clear distinction between batches. I'll update the question with the figure. Can you please elaborate more on calculating a factor to be taken into account for later steps? Do you mean to extract the principal component correlated with batch and later do something with it?

Well, what I do is include the known batch effects on the linear models. That could be via the component of the PCA or the categories you know of your batches. Could you expand into what are your batches ? What are 3_8, 3_38 ... (I assume M is male and F is Female).?

@Llopis unfortunately I have no metadata information regarding the first part of the mouse id. They only provide sex (M/F) which is as you assumed.

Your PC2 separates cells of the mice but is only 0.6% of the variation, so I would say there isn't batch effect. The first dimension is quite high, but I don't know if this is normal in scRNA-seq. I wouldn't adjust nor remove batch effect here if this would be RNA-seq. But I have never analysed scRNA

Let us continue this discussion in chat.

in my experience, the absolute first thing you need to do is normalize for library size. I suspect that if you color your cells according to size you will notice a clear correlation with PC1.

MAGIC assumes input data has been both library-size normalized, and either log or sqrt transformed prior to imputation (see also: MAGIC tutorial). Additionally, any graph-based methods (MAGIC, PHATE, t-SNE, UMAP, spectral clustering, Louvain, etc etc) will give flawed results if your data contains a batch effect, since the neighbourhood graph would reflect that structure of your batch effect, and worse, imputation would further reinforce this batch effect.
Thus I would recommend the following pipeline:

Library-size normalization
Square root (or log) transform
Batch effect removal
Imputation

Regarding your update, the reason you don't see the batch effect in the raw counts data is simply that the batch effect is not visible in the most highly expressed genes. Prior to transformation, the principal source of variation in your data is simply the expression of the most highly expressed genes - this is essentially masking the source of the batch effect, not removing it. I recommend never working with raw molecule counts in scRNAseq, as the raw counts data hides much of the heterogenity in your dataset, which is precisely what you are looking in doing single-cell RNA-seq.


text:
Is there anything similar to GSEA for locus-based (instead of of gene-based) data?
As the question states, I am interested in an analysis similar to Gene Set Enrichment Analysis (ranked gene sets) but focused on locus-level data instead of genes.
To explain in greater detail: I have a set of genomic coordinates from DNA methylation data that are scored based on their contribution to components from a non-negative matrix factorization. I am interested in better understanding the biological underpinnings of what these NMF components represent.  I know I can use Locus Overlap Analysis tools like GREAT or LOLA to tell me what GO/Functions are enriched in a particular BED file, but this type of analysis doesn't take into account the scores/rank of the loci in the BED file as would be done in an analysis akin to GSEA.
Does anyone know of any tools/methods for analyzing functional enrichment while taking advantage of feature/locus ranks/scores with locus-level data?
I know I probably didn't explain this very clearly, so please let me know if I can clarify anything for you.  Thanks

I tinkered with a program GSEA-SNP quite a few years back, which claims that it does a similar ranking procedure with SNPs. It carries out its procedure by first linking SNPs with genes, then running an algorithm similar to GSEA. There's a bit more detail in the paper.
Unfortunately this is a space I've got a bit further away from in recent years, so I don't know if anyone has carried out the legwork to generalise the algorithm for sets of any type where each item within each set fits a continuous (but not necessarily normal) distribution.

This is a step in the right direction and better than I have been able to find, so thank you!  I haven't been able to find a non GWAS version of these tools, but these versions might work for my application. Thanks!

You can use GSEA-like methods for any type of data. The basis of the test is to analyze if there is a preference for a group (a gene set) in the ranking of the other variable.
So you first need to find where the functions you are interested are in the genome (To associate a range/position with a given function). Then using this association you can test if the ranking of the locus is associated with a specific function that is in a given area.
I haven't seen anything like that implemented but you could do it yourself using biomaRt to retrieve data from where a gene (or a function) is located and fgsea to test which functions are associated with your score.


text:
determining doublets in single-cell RNA-seq
Doublets are a known problem with scRNA-seq experiments, where 2 or more cells are sometimes captured instead. To determine their presence, there are studies that mix multiple species (such as human and mouse) or distinct cell types and then count how many cells contain genes from the very different transcriptomes. In those experiments, you can pick any cell and determine if it's a doublet based on a few key genes. In a typical single-cell experiment on one sample, you don't have the luxury of knowing what cells to expect and if a particular cell looks abnormal.
If you have a large fraction of doublets, you should be able to see a double peak if you plot the number of genes/transcripts per cell. However, if you have something like 1% doublets, you probably wouldn't notice anything. It seems that most people simply filter out top X% cells based on the number of genes/transcripts. How do you know what X should be? Can you actually estimate the fraction of doublets in a given experiment?
Update:
I found a nice figure from Stoeckius et al. to illustrate the issue:

Although the number of transcripts/UMIs seems to be roughly double for doublets/multiplets (median of ~1000 vs ~500), there is still high overlap.

related question here, which asks what doublets are, rather than how to identify them.

actually, that question inspired mine

Expected rates of doublets / duplets / multiplets

Fluidigm C1 doublet rate: around 1-5% depending on chip type used.
More information: Fluidigm white paper: Redesign of C1 Medium-Cell 96 IFCs Improves Single-Cell Capture Efficiency (PN 101-3328 B1)
10x Genomics Chromium 3'mRNA-Seq doublet rate: around 1% per 1000 cells captured

Rates in the 'Chromium Single Cell 3’ Reagent Kits v2 User Guide':
Multiplet Rate (%)  # of Cells Loaded   # of Cells Recovered
0.40%   870 500
0.80%   1700    1000
1.60%   3500    2000
2.30%   5300    3000
3.10%   7000    4000
3.90%   8700    5000
4.60%   10500   6000
5.40%   12200   7000
6.10%   14000   8000
6.90%   15700   9000
7.60%   17400   10000

Detection and identification of doublets and multiplets
Experimental

Visual inspection (e.g. Fluidigm C1)
FACS index sorting data can be used to identify doublets.
Mixed species experiment.

Illumina | BioRad White paper Single-Cell RNA Data Analysis Workflow
Zheng et al., 2017: Massively parallel digital transcriptional profiling of single cells

Cell "hashing" with oligo-tagged antibodies (Stoeckius et al., 2017): cells (barcodes) positive for two or more oligos are multiplets.

Computational

Mixed gene expression (e.g. lymphoid and myeloid, or T and B cell)
Ilicic et al. describe an approach for detecting low quality cells (incl. multiplets). They train a support vector machine (SVM) model on a dataset annotated by microscopy inspection. They find, among others, that

"wells containing multiple cells (multiples) show similar expression and noise patterns to broken cells", suggesting "that multiple cells contain a mixture of broken and high quality cells."
difference in proportion of duplicated reads
higher ERCC spike-in ratios

demuxlet (Kang et al., 2017, Nature Biotechnology)
DoubletDetection

DoubletDetection looks very promising. It seems to be the only method that does not require any prior knowledge about the cells or a mixture of distinct cell types/samples.

We cannot assume that doublets will produce more UMIs
I would caution the assumption that all doublets will have twice the UMI levels of isolated single-cells. Many "doublets" could contain multiplets of 3 or more cells, depending on how many cells have been loaded into the experiment. Most vendors of single-cell technologies quote multiplet rates for doublets, triplets, etc combined.
Furthermore, we cannot assume that the data from multiplet (doublet) cells is of high quality or produces comparable reads and UMI counts than we would expect from two single cells combined. Some tests of doublets (such as using mixed mouse and human samples) have found that the data generated from doublet droplets or nanowells has relatively low UMI abundance. This may be due to the protocols having been optimised for droplets or nanowells containing only one cell. Especially where tissue preparation is required this is a concern and is consistent with observations that doublet cells appear to be broken or contaminated with debris when they can be identified. This is highly dependent on the tissue you and working with, sample preparation, and the platform used to separate single-cells.
I have not been able to reproduce the above result using data from various species (using total UMI counts per droplet or housekeeper genes), although this is work still in progress. Even in the data shown, many multiplet samples have similar total UMI counts to the true single-cells. Therefore, they cannot be separated by a simple threshold (despite being significantly different). There have been several tools developed already for detecting multiplets (and filtering them out) but this is still an area of active development and there are more ways that these could be refined. In many single-cell analyses, arbitrary thresholds are being used to filter the data before further analysis and this may be inadvertently reducing the impact of multiplets if they have produced fewer UMI counts than single-cells for which the protocols have been optimised for.
Experimental protocols can be optimised to reduce doublets
If you have any insight into how these bioinformatics solutions to issues with multiplets could be used or improved I encourage you to pursue it. However, I suspect that the eventual solution to this issue is in optimising the single-cell NGS experiments for the sample type and platform that you are working with. Technology in single-cell 'omics is continuing to develop and any solution which produces more useful single-cell data than contaminants or doublets (that need to be removed), the more cells can be sampled at the same cost. There is a strong incentive for experimentalists and single-cell technology vendors to optimise this issue.


text:
In GFF3, annotating more than one protein-coding gene (i.e. polycistronic) contained in a eukaryotic mRNA
Below is an example of a simple GFF3 file:
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS 3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS 3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

My question is: if we found another coding sequence (encoding a different protein) range from 3752 to 3904, how should the GFF3 be formatted? It seems to me that the GFF3 file can only allow one protein-coding gene per mRNA. If not, could anyone show me one example? Thank you!
[Edit]:
Sorry, using "CDS" is misleading. I changed the title and contents to "contain more than one protein-coding genes" instead of "more than one CDS". I also change the starting position of the second protein to 3752 (was 3750) to reflect the fact. My sincere apologies.

The GFF3 file format specification doesn't care where annotations are described in the file, although it can make some things easier if the annotations are ordered by the start position:
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3750    3904    .   -   0   ID=ATRV01011:CDS:1;Parent=ATRV01011.0
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

There are many examples in well-studied genomes of overlapping genes, so the GFF3 format needs to be flexible enough to handle these. Different genes can be encoded across the same region, but on different strands. But different genes can also be encoded on the same strand, sometimes sharing an exon.

thanks for your help. So essentially we can just use a new name (e.g here from ATNG to ATRV) to describe the new protein, correct? Also please see my edit. The second protein should actually be in a different open reading frame. Sorry I calculate it wrong the first time. However, I assume your answer still apply.

It is certainly possible to encode polycistrons in GFF3. In fact, this is discussed explicitly in the GFF3 specification.

Polycistronic transcripts
This is the case in which a single (possibly spliced) transcript encodes multiple open reading frames that generate independent protein products.
----->XXXXXXX*-->BBBBBB*--->ZZZZ*-->AAAAAA*-----
Since the single transcript corresponds to multiple genes that can be identified by genetic analysis, the recommended solution here is to create four "gene" objects and make them the parent for a single transcript. The transcript will contain a single exon (in the unspliced case) and four separate CDSs:
chrX  . gene XXXX YYYY  .  +  . ID=gene01;name=resA
chrX  . gene XXXX YYYY  .  +  . ID=gene02;name=resB
chrX  . gene XXXX YYYY  .  +  . ID=gene03;name=resX
chrX  . gene XXXX YYYY  .  +  . ID=gene04;name=resZ
chrX  . mRNA XXXX YYYY  .  +  . ID=tran01;Parent=gene01,gene02,gene03,gene04
chrX  . exon XXXX YYYY  .  +  . ID=exon00001;Parent=tran01
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene01
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene02
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene03
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene04

To disambiguate the relationship between which genes encode which CDSs, you may use the Derives_from relationship.

This example doesn't indicate what the coordinates of each feature should be, but I think it makes most sense to have the gene and CDS feature coordinates correspond to the region coding each protein product, and the mRNA and exon feature coordinates correspond to the entire span.
BIG CAVEAT: Polycistrons are not typically considered a normal feature of eukaryotic genomes—they're much more common in bacteria and archaea. However, there are some confirmed cases of eukaryotic polycistronic mRNAs. That said, if overlapping gene predictions are observed in a GFF3 file, one cannot immediately assume that these are polycistronic. One or both of the gene predictions might be technical artifacts, or real genes that are incorrectly annotated, or (more rarely) actual overlapping genes. A claim that multiple protein products are expressed from a single transcript would require some convincing experimental validation.

If you found another CDS from 3750 to 3904 then it would necessarily be on a different transcript, given that you already have a CDS from 3860 to 3913 and having overlapping CDSs on the same transcript doesn't make much biological sense (though I'm sure there's some weird organism out there that has a gene like this). So the likely GFF3 file would look like:
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.2;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010.2:exon:1;Parent=ATNG01010.2
1   T1  CDS     3759    3904    .   +   0   ID=ATNG01010.2:CDS:1;Parent=ATNG01010.2
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010.2:exon:2;Parent=ATNG01010.2
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010.2:CDS:2;Parent=ATNG01010.2
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010.2:exon:3;Parent=ATNG01010.2

Note that I've just added another transcript, since that's likely what would be biologically going on.

I calculate the ORF wrong (should be from 3752 to 3904, not 3750 to 3904). Please see the updated question. It does overlap with the original CDS but not in the same frame. This is a type of uORF. Any thought? Again, thanks for your help and sorry for the mistake.

Then the answer from @gringer is the more relevant one.

I found one similar example in an Arabidopsis annotation file. It is case of the cpuORF (CONSERVED PEPTIDE UPSTREAM OPEN READING FRAME ). Although the first and the third mRNA isoforms are identical, they repeated the mRNA isoform and name them differently to accommodate the additional ORF/CDS. I guess this is another way to do it.

See the ranges for AT5G03190.1 isoform and AT5G03190.3 isoform.
Chr5    TAIR10  gene    758374  760382  .   +   .   ID=AT5G03190;Note=protein_coding_gene;Name=AT5G03190
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.1;Parent=AT5G03190;Name=AT5G03190.1;Index=1
Chr5    TAIR10  protein 758793  760148  .   +   .   ID=AT5G03190.1-Protein;Name=AT5G03190.1;Derives_from=AT5G03190.1
Chr5    TAIR10  exon    758374  760382  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  five_prime_UTR  758374  758792  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  CDS 758793  760148  .   +   0   Parent=AT5G03190.1,AT5G03190.1-Protein;
Chr5    TAIR10  three_prime_UTR 760149  760382  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.2;Parent=AT5G03190;Name=AT5G03190.2;Index=1
Chr5    TAIR10  protein 758539  760148  .   +   .   ID=AT5G03190.2-Protein;Name=AT5G03190.2;Derives_from=AT5G03190.2
Chr5    TAIR10  exon    758374  758660  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  five_prime_UTR  758374  758538  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  CDS 758539  758660  .   +   0   Parent=AT5G03190.2,AT5G03190.2-Protein;
Chr5    TAIR10  exon    758843  760382  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  CDS 758843  760148  .   +   1   Parent=AT5G03190.2,AT5G03190.2-Protein;
Chr5    TAIR10  three_prime_UTR 760149  760382  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.3;Parent=AT5G03190;Name=AT5G03190.3;Index=1
Chr5    TAIR10  protein 758539  758676  .   +   .   ID=AT5G03190.3-Protein;Name=AT5G03190.3;Derives_from=AT5G03190.3
Chr5    TAIR10  exon    758374  760382  .   +   .   Parent=AT5G03190.3
Chr5    TAIR10  five_prime_UTR  758374  758538  .   +   .   Parent=AT5G03190.3
Chr5    TAIR10  CDS 758539  758676  .   +   0   Parent=AT5G03190.3,AT5G03190.3-Protein;
Chr5    TAIR10  three_prime_UTR 758677  760382  .   +   .   Parent=AT5G03190.3


text:
ROC analysis, CAFA 2 experiment
My question is about ROC curves used in CAFA2 experiment. In this paper they used the ROC analysis for the term-centric evaluation. In order to
perform ROC curve analysis we should have a continuous variable and a classifier (categorical) variable. In CAFA2 experiment,
which variables did they consider as a categorical and continuous? Thank you.

The paper doesn't seem to have been written for general readers. Unfortunately I don't understand the equations enough to be able to explain what the categorical and continuous variables are. Just in case it helps, here are the equations from the paper showing the precision/recall calculations:
$$\begin{array}{@{}rcl@{}} \text{pr}(\tau) &=& \frac{1}{m(\tau)}\sum\limits_{i=1}^{m(\tau)} \frac{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in P_{i}(\tau) \wedge f \in T_{i}\right)}{\sum_{f} \unicode{x1D7D9}\left(f \in P_{i}(\tau) \right)},\\ \text{rc}(\tau) &=& \frac{1}{n_{e}}\sum\limits_{i=1}^{n_{e}} \frac{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in {P}_{i}(\tau) \wedge f \in T_{i}\right)}{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in T_{i} \right)}, \\ F_{\max} &=& \max_{\tau} \left\{ \frac{2\cdot \text{pr}(\tau)\cdot \text{rc}(\tau)}{\text{pr}(\tau) + \text{rc}(\tau)} \right\}, \end{array} $$

where P i (τ) denotes the set of terms that have predicted scores greater than or equal to τ for a protein sequence i, T i denotes the corresponding ground-truth set of terms for that sequence, m(τ) is the number of sequences with at least one predicted score greater than or equal to τ, $\unicode{x1D7D9}\left (\cdot \right)$ is an indicator function, and n e is the number of targets used in a particular mode of evaluation.

From a quick reading of the paper it looks like they are predicting biological function of proteins. They frame this as a classification problem, where they try to classify a protein as pertaining to a set of terms (Abnormality of the lymph nodes, Abnormality of skin physiology, etc). It looks like they predict for each protein whether it falls into this category and create a ROC curve to plot the false positive and true positive rates of prediction.
In this case, the continuous variable would be your predictions (Abnormality of the lymph nodes: predicted 0.65 probability, Abnormality of skin physiology: predicted 0.25 probability). The categorical variable would then be the ground truth that you are trying to predict (Abnormality of the lymph nodes: known 1, Abnormality of skin physiology: known 0). Based on these two variables you can make your ROC curve.


text:
Why doesn't Biopython AlignIO.read() recognise the 'mauve' format?
On this Biopython tutorial, they describe how to import a multiple sequence alignment in the Mauve (XMFA: extensible multi fasta format). So I imported the AlignIO module:
from Bio import AlignIO
alignment = AlignIO.read(open("alignment.xmfa"), "mauve")

but I get the following error:
ValueError
Traceback (most recent call last) <ipython-input-30-396e54791d7a> in <module>()
----> 1 alignment = AlignIO.read(open("/Users/cr517/Documents/alignment.xmfa"),
"mauve")

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in read(handle, format, seq_count, alphabet)
    422     iterator = parse(handle, format, seq_count, alphabet)
    423     try:
--> 424         first = next(iterator)
    425     except StopIteration:
    426         first = None

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in parse(handle, format, seq_count, alphabet)
    366                                                 seq_count=seq_count)
    367         else:
--> 368             raise ValueError("Unknown format '%s'" % format)
    369
    370         # This imposes some overhead... wait until we drop Python 2.4 to fix it

ValueError: Unknown format 'mauve'

You're using version 1.68 or older. Mauve support was added in 1.70.


text:
How do I compare transcriptomes of the same species across different developmental stages?
I want to identify the most highly expressed genes across several developmental stages of a given insect, in hopes of identifying a protein that can be used to design an insecticide.
I have a few doubts :
Objective

Several studies have been performed on individual developmental stages and their transcriptome has been published. Should I compare each set of data against the other manually, to find genes that are expressed universally and then target the most highly expressed ones or is there any software / computational method that can be used instead?

Question
Is it right to compare transcriptomes to find commonly expressed genes?
I did try to look at the materials and methods section of a few papers on comparative transcriptomics, but it only says that data was compared without stating the method. One paper cited VennBLAST, but I wasn't able to find the tool online.

Could you clarify what do you mean by individual developmental stages? They looked at the same individual in several stages or they look in several samples different stages of development. What I think you want to do is find differentially expressed genes, search for those terms and let us know if you need something else

I mean that each paper concentrates on a single developmental stage. They analyze several samples for different stages of development.
Also, no, I do not want to find differentially expressed genes.From what I understand, differentially expressed genes are those that are upregulated only in a particular cell type (in my case stage). I want to find out if a common set of genes are being upregulated throughout the life cycle of the insect, ie. they are common to all the life stages. I thought I can do this by analyzing the transcriptomes of each stage. I hope this is clear!

You can find the differentially expressed genes that are downregulated too, and it is not specific to a cell type that depends on your data. But it is a comparison between two stages, so you can compare all the initial stage with all the following stages. Then you find the common genes that are up-regulated in all the cases of interest

Oh! I didn't know that differentially downregulated genes can be found too! Thanks for your answer :D

This was a comment!! You seem to need more guidance, ask your advisor or a senior people how to find the most expressed genes and from there how to find it those that are so in several stages.