text:
Fitting branch models with godon
Is it possible to fit branch models (Yang and Nielsen, 1998) in godon? That is, the model which averages over all sites, as opposed to the branch site models?

For those VTCing this - Bioinformatics.SE is the specified help forum for the godon package. If you've got a problem, take it up with the meta

The votes to close are because your question lacks a lot of detail.

And yet, it has all of the information necessary for the question to be answered.

Voted to open, because it's exactly the type of questions we aim for - concise, concrete and answerable.

Branch-models are currently not implemented in godon. In the PAML manual Ziheng Yang writes:

This model is very parameter-rich and its use is discouraged.

This was one of the reasons I did not implement it.
Disclaimer: I am the author of godon.


text:
Biopython Phylogenetic Tree replace branch tip labels by sequence logos
Having recently constructed a lot of phylogenetic trees with the module TreeConstruction from Phylo package from Biopython, I've been asked to replace the branch tip labels by the corresponding sequence logos (which I have in the same folder). I thought that it would be more efficient to make a code to generate the logo-trees automatically, as I would have to make a lot of them.
The first idea that I came up with was to see whether the functions used to build the tree had an argument to replace the branch tip labels or to remove them, which I couldn't find. Therefore, I removed the branch tip labels by setting their font size to 0: (following is the code to build the tree)
# Modules to build the tree
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo import draw
from Bio import Phylo, AlignIO
import subprocess
import matplotlib
import matplotlib.pyplot as plt

alignment = AlignIO.read('MotifSeqAligned.fasta', 'fasta') # reading the alignment file

calculator = DistanceCalculator('ident')
dm = calculator.get_distance(alignment) # distance matrix

constructor = DistanceTreeConstructor()
tree = constructor.nj(dm) # build with neighbour joining algorithm a tree from dm

Phylo.write(tree, 'TreeToCutOff.nwk', 'newick')

plt.rc('font', size=0)          # controls default text sizes #HERE IS THE SETTING FOR THAT ALLOWS ME TO HIDE THE BRANCH TIP LABELS
plt.rc('axes', titlesize=14)     # fontsize of the axes title
plt.rc('xtick', labelsize=10)    # fontsize of the tick labels
plt.rc('ytick', labelsize=10)    # fontsize of the tick labels
plt.rc('figure', titlesize=18)   # fontsize of the figure title

draw(tree, do_show=False)
plt.savefig("TreeToCutOff.svg", format='svg', dpi=1200)

From this code I could get the tree:

Since I don't know how to get the y coordinates of the branches to add the logos one by one, I built a column of logos with matplotlib, that I intended then to paste on the tree in python. The code to build the column of logos is the following:
#Extract filename from newick
newickFile = open("TreeToCutOff.nwk", 'r').read()
orderedLogos = ["{}.eps".format(i) for i in re.split('(\W)', newickFile) if "Profile" in i]

#Initialize the figure
fig = plt.figure()

#Add each image one after the other in the right order
for i, files in enumerate(orderedLogos):
    img1 = mpimg.imread(files)
    ax1 = fig.add_subplot(len(orderedLogos), 1, 1+i)
    ax1.imshow(img1)
    ax1.set_xticks([])
    ax1.set_yticks([])

# plt.show()
plt.savefig("RowsOfLogos.svg", format='svg', dpi=1200)
plt.clf()
plt.cla()

Having my tree and the column of logos in .svg or .png, I couldn't find any way to stack them properly. My first idea was to use the library svgutils which seemed to be easy to handle, with the following code: (taken from svgutils tutorials)
import svgutils.transform as sg
# Assemble
#create new SVG figure
fig = sg.SVGFigure("14cm", "14cm")

# load figures
fig1 = sg.fromfile('TreeToCutOff.svg')
fig2 = sg.fromfile('RowsOfLogos.svg')

# get the plot objects
plot1 = fig1.getroot()
plot2 = fig2.getroot()
plot2.moveto(280, 100, scale=0.05)

# append plots and labels to figure
fig.append([plot1, plot2])

But the issue with the output was that the background of the column of logos was white and thus, I was pasting a huge white image with a thin column of logos on the tree. And I couldn't find a way to crop the column of logos with svgutils. I tried the module Image from PIL package to build a tree of logos from .png files but couldn't see the tree used as background after pasting the column of logos.
There may be a way to do what I'm aiming for with matplotlib (which would be to stack 2 files .png, and place the logos all the time at the same distance), but I couldn't work it out.
Does anyone know what the best solution is to make a tree of logos (as in the following image which I could only build manually with inkscape) with python libraries, allowing to automate the process without having to adapt the code depending on the number of branches?

Following is a subset of "MotifSeqAligned.fasta" containing the aligned sequences used to build the trees:
>ProfileCluster0.meme
---SSNDTTTCCAGGAAD-
>ProfileCluster1.meme
YBNRD---TTCYYGGAAT-
>ProfileCluster10.meme
---VDKDWTTCTYGGAAT-

With the 3 corresponding logos:

The full "MotifSeqAligned.fasta" and all the logos.eps (as I used them instead of .png which is the format asked by the forum) can be found here.

Thank you for putting so much effort into your question! If I could ask for one more thing though, it would be great if you could give us a minimal, complete and verifiable example, code that we can use directly on our system to replicate your issue so we can play with it and see if we can find a solution.

Thank you for the suggestion, I added the inputs at the end of the topic. Note that the logos are .eps files and here I added them as images (.png). The full inputs are on the link at the end of the topic.

If there's no obvious biopython solution have a look at the ETE Toolkit

To expand on my comment from yesterday. You could do this with the ETE Toolkit (I just copied one logo file rather than converting all 26 to png):
from ete3 import Tree, TreeStyle, faces
def mylayout(node):
    if node.is_leaf():
        logo_face = faces.ImgFace(str.split(node.name, '.')[0] + ".png") # this doesn't seem to work with eps files. You could try other formats
        faces.add_face_to_node(logo_face, node, column=0)
    node.img_style["size"] = 0  # remove blue dots from nodes

t = Tree("tree.nwk", format=3)
ts = TreeStyle()
ts.layout_fn = mylayout
ts.show_leaf_name=False # remove sequence labels
ts.scale =  10000 # rescale branch lengths so they are longer than the width of the logos
t.render("formatted.png", tree_style = ts, h=3000, w=3000) # you may need to fiddle with dimensions and scaling to get the look you want

If you want all of the logos lined up in a column add aligned=True to faces.add_face_to_node

Thanks! That's exactly what I was looking for! I will try with all my logos and let you know if I have further issues.

source: Stack Exchange

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

sort: votes


text:
Retrieve ID ligand from PDB file
I have thousands of PDB files and I need to extract the following information from each of them: an ID or the name of the small molecule (ligand) present in the complex.
Is there a way to do that with BioPython or similar tools? Taking this PDB entry as an example, the tools should return 7FM.

Can you please [edit] your question and include some examples and the output you would want from them? That way, we can test our solutions using your examples and we can be sure to give you the output you need.

The PDB entry states that the ligand is 7FM.

Anyway use the MMCIF2Dict tool to extract  the information from the header via the PDB module of BioPython. If it's not in the header you can't do it because XML is not supported (not currently supported). Goodluck

And how do you use this function?

question is not clear anymore is the author referring to files a s PDB format or PDBx/mmcif format ? In case of the former Biopython can get the name of an HETATM residue see https://files.rcsb.org/view/5TO8.pdb

I found out that the package PyPDB has the function get_ligands that does exactly what I'm looking for.

Thanks, for the feedback. Apologies for not spotting this.

No problem! Strange I was not able to find this tool before...

copying from How to download ligands for PDB structure one of the answers, one of mines, code  :

import requests

import json

parsed_lig_dict = {}

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    # print(risposta.status_code)
    # print((parsed))
    # print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    # print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    # print('\n\n',parsed_lig, type(parsed_lig))

    # for i in parsed:
    #     print('\n', i  , parsed[i])

    for lig in parsed_lig:
        risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity/{entry_id}/{lig}")
        parsed = json.loads(risposta.content.decode())
        # print('\n',risposta.status_code, risposta.content)

        # print(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"],
              # parsed["pdbx_entity_nonpoly"]["comp_id"])

        parsed_lig_dict[lig] = [parsed["pdbx_entity_nonpoly"]["comp_id"] ,
                                dict.fromkeys(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"])]

    for lig in parsed_lig_dict:
        for chain in parsed_lig_dict[lig][1]:
            risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity_instance/{entry_id}/{chain}")
            # print(risposta.status_code)

            parsed = json.loads(risposta.content.decode())
            # print('\n',risposta.status_code, risposta.content)
            # print(chain , '--> ',parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"])

            parsed_lig_dict[lig][1][chain] = parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"]

    # print(parsed_lig_dict)

    # for i in parsed_lig_dict:
    #     print(parsed_lig_dict[i])

    # for i in parsed_lig_dict:
    #     for y in parsed_lig_dict[i][1]:
    #         print(i,parsed_lig_dict[i][0], y , parsed_lig_dict[i][1][y])

    cnt = 1
    for i in parsed_lig_dict:
        for y in parsed_lig_dict[i][1]:
            seq_id = parsed_lig_dict[i][1][y]
            chain = y
            comp_id = parsed_lig_dict[i][0]
            print(seq_id , y)

            risposta = requests.get(f"https://models.rcsb.org/v1/{pdb_file}/ligand?auth_seq_id={seq_id}&label_asym_id={chain}&encoding=sdf" ,
                              allow_redirects=True)
            print(risposta.status_code)

            filename = f"lig_{comp_id}_{chain}_{seq_id}_{cnt}.sdf"
            open(filename, 'wb').write(risposta.content)
            cnt += 1

if __name__ == '__main__':

   pdb_file = '5to8'

   main(pdb_file)

output file : lig_7FM_B_701_1.sdf :
7FM
  ModelServer 0.9.10

 39 43  0  0  1  0  0  0  0  0  0
   -2.0610  -10.3590   40.5700 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.1170  -11.7950   42.6720 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8640  -12.6100   43.6800 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6970  -10.3290   36.5750 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8210  -10.9080   35.9890 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6600   -8.1940   35.9530 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.0300   -9.0340   41.4650 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.9390   -9.6680   41.0120 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7980  -10.0800   39.7250 N   0  0  0  0  0  0  0  0  0  0  0  0
    1.7720   -9.8430   38.8300 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.9200   -9.1930   39.2330 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.9910   -8.8020   40.5580 C   0  0  0  0  0  0  0  0  0  0  0  0
    4.0340   -8.8690   38.3150 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0830   -9.8810   41.9510 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3380  -10.4670   41.7760 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.3020  -10.9670   40.4240 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8290  -11.6870   41.4760 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8800  -11.1950   42.8370 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.0400  -12.4130   41.5840 N   0  0  0  0  0  0  0  0  0  0  0  0
   -5.2460  -12.6990   43.0180 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.5160  -10.2930   37.5490 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.4060  -11.1430   37.2090 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.8420  -10.2200   35.4050 N   0  0  0  0  0  0  0  0  0  0  0  0
   -2.7360   -8.8690   35.4000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6240   -8.9330   36.5380 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.9670  -12.3410   35.9950 N   0  0  0  0  0  0  0  0  0  0  0  0
   -2.5680  -13.0870   37.1550 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.4690  -13.1930   34.5670 S   0  0  0  0  0  0  0  0  0  0  0  0
   -4.0440  -12.7310   37.3480 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.9670  -12.8650   40.6430 C   0  0  0  0  0  0  0  0  0  0  0  0
   -6.9720  -13.5050   40.9690 O   0  0  0  0  0  0  0  0  0  0  0  0
   -4.6200  -13.4130   38.5900 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.8010  -12.6200   39.1510 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.5520  -12.3690   33.7940 O   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1430  -14.5640   34.8900 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.0530  -13.2140   33.7590 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.0490   -8.1440   38.8710 F   0  0  0  0  0  0  0  0  0  0  0  0
    4.6120   -9.9430   37.7390 F   0  0  0  0  0  0  0  0  0  0  0  0
    3.5350   -8.1080   37.3100 F   0  0  0  0  0  0  0  0  0  0  0  0
  1 16  2  0  0  0  0
  1 15  1  0  0  0  0
  2 17  2  0  0  0  0
  2 18  1  0  0  0  0
  2  3  1  0  0  0  0
  3 20  1  0  0  0  0
  4  5  2  0  0  0  0
  4 25  1  0  0  0  0
  4 22  1  0  0  0  0
  5 23  1  0  0  0  0
  5 26  1  0  0  0  0
  6 24  1  0  0  0  0
  6 25  2  0  0  0  0
  7 12  2  0  0  0  0
  7  8  1  0  0  0  0
  8  9  2  0  0  0  0
  8 14  1  0  0  0  0
  9 10  1  0  0  0  0
 10 21  1  0  0  0  0
 10 11  2  0  0  0  0
 11 13  1  0  0  0  0
 11 12  1  0  0  0  0
 13 39  1  0  0  0  0
 13 38  1  0  0  0  0
 13 37  1  0  0  0  0
 14 15  1  0  0  0  0
 15 18  2  0  0  0  0
 16 17  1  0  0  0  0
 17 19  1  0  0  0  0
 19 30  1  0  0  0  0
 19 20  1  0  0  0  0
 21 22  1  0  0  0  0
 23 24  2  0  0  0  0
 26 28  1  0  0  0  0
 26 27  1  0  0  0  0
 27 29  1  0  0  0  0
 28 36  1  0  0  0  0
 28 34  2  0  0  0  0
 28 35  2  0  0  0  0
 29 32  1  0  0  0  0
 30 33  1  0  0  0  0
 30 31  2  0  0  0  0
 32 33  1  0  0  0  0
M  END
> <model_server_result.job_id>
RFDalM1AmREvgkv0PqdfQw

> <model_server_result.datetime_utc>
2023-10-27 22:13:14

> <model_server_result.server_version>
0.9.10

> <model_server_result.query_name>
ligand

> <model_server_result.source_id>
pdb-bcif

> <model_server_result.entry_id>
5to8

> <model_server_params.name>
atom_site

> <model_server_params.value>
{"label_asym_id":"B","auth_seq_id":701}

> <model_server_stats.io_time_ms>
5

> <model_server_stats.parse_time_ms>
76

> <model_server_stats.create_model_time_ms>
2

> <model_server_stats.query_time_ms>
244

> <model_server_stats.encode_time_ms>
0

> <model_server_stats.element_count>
39

$$$$

picture of .sdf file :

my attempt, will get only non water (H2O) or non aa (aminoacid) present in the PDB .... tricky if the ligand is an aa :
from Bio.PDB import PDBParser

from Bio.PDB import Selection

d = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
     'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
     'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
     'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

pdb_file = '5to8'

struct = PDBParser(QUIET = True).get_structure(pdb_file, pdb_file+'.pdb')

res_list = Selection.unfold_entities(struct , 'R')

for i in res_list :

    if i.resname != 'HOH' and i.resname not in d :

        print(i, i.resname)

output :
<Residue 7FM het=H_7FM resseq=701 icode= > 7FM

authors: Юлия Иванькова

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


text:
Filtering out reads from a reference (e.g. rRNA) using bowtie2
As I understand it, bowtie2 can easily be used to split reads into one of two groups:

reads for which both of a pair align well to a reference (using e.g. --al-conc-gz)
reads for which one or both of a pair do not align a reference (using e.g --un-conc-gz)

But I really want to split this second group into reads for which neither of a pair align to the reference.
bowtie2's report give a nice breakdown, but it's not clear to me how to extract specific members of these groups:
1000000 reads; of these:
  1000000 (100.00%) were paired; of these:
    958118 (95.81%) aligned concordantly 0 times
    1329 (0.13%) aligned concordantly exactly 1 time
    40553 (4.06%) aligned concordantly >1 times
    ----
    958118 pairs aligned concordantly 0 times; of these:
      179 (0.02%) aligned discordantly 1 time
    ----
    957939 pairs aligned 0 times concordantly or discordantly; of these:
      1915878 mates make up the pairs; of these:
        711518 (37.14%) aligned 0 times
        32134 (1.68%) aligned exactly 1 time
        1172226 (61.18%) aligned >1 times
64.42% overall alignment rate

In this case, I really want those that didn't align at all (i.e. the 35.58% not included in the 64.42% overall alignment rate which is the same as the 37.14% of the pairs that did not align concordantly).
Since in my case the reference is rRNA and many of those sequences are shorter than my reads, I'm not surprised that so there are so few concordant reads. I'm really interested in reads that aren't concordant and aren't discordant discordant either.
So, if I'm sticking to bowtie2, it looks like my only option is to take the resulting SAM/BAM file and filter what I want based on flags and then convert that back into FASTQ files.
Is there another option that I'm missing? It seems that the remaining output options relate to unpaired reads or SAM files.
I'm looking into other aligners now...
My goal is to go from FASTQ files to FASTQ files without extra steps, but I can go through a SAM/BAM file stage if I need to and even build a solution around that if needed.

If you want to remove rRNA, use tools that are made for this. We usually go for SortMeRNA, which can do exactly what you describe your goal is, i.e. go from FASTQ to FASTQ, splitting your reads into rRNA and non-rRNA.
In our experience, it is very accurate, if a bit slow at times.

Bam flags will indicate which reads have both pairs not mapping.

Nothing mapped.  You need to figure out why, and you can look at just about anything out figure that out.

Yes, thanks for mentioning SAM/BAM flags. I thought I mentioned them, but it looks like I only implied them. (Now updated to mention them)

I'm trying to filter out rRNA reads. So about 2/3 of the read pairs had at least one read align to my rRNA reference. However, I want the remaining 1/3 of the read pairs for which neither read aligned at all.

You are trying to filter out rRNA read by aligning them to a reference of rRNAs?  This seems ill-advised.  Filter out rRNA after aligning to your whole genome.


text:
Alignment with inserts and keeping the indexing of ref seq intact
Parts of sequences are given below-
Reference sequence (pre-alignment):
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAAAATCTGTGT

Reference sequence (post-alignment) and below it is Sample sequence (post-alignment):
--------------------------------------------------------------------------------------------------attaaa---------ggtt------------------tataccttc---------ccaggtaacaaa-------------ccaacc-----aactttcgatctcttgtagatctgttctctaaacgaactttaaaatctgtgt
--------------------------------------------------------------------------------------------------------------------------------------------------------------taacaaa-------------ccaacc-----aactttcgatctcttgtagatctgttctctaaacgaactttaaaatctgtgt

I'm adding a simpler to interpret example-
Say the ref seq is aattaaatttgggggtttt and the sample seq is ttaaggggttaaatttgggggt--t. Then post-alignment, they will be like-
--aa----ttaaatttgggggtttt
ttaaggggttaaatttgggggt--t

Now, since my ref seq was
0    1    2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18
a    a    t    t    a    a    a    t    t    t    g    g    g    g    g    t    t    t    t

I want that post-alignment also, the indexing should be conserved-
           0    1                        2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18
 -    -    a    a    -    -    -    -    t    t    a    a    a    t    t    t    g    g    g    g    g    t    t    t    t
 t    t    a    a    g    g    g    g    t    t    a    a    a    t    t    t    g    g    g    g    g    t    -    -    t

I want to keep the indexing of the reference sequence conserved(i.e. the first base in ref seq post-alignment is a, second is t, third is t, etc.), like they do in standard softwares, and then I want to run some quick analysis on it, say to check for conservation of a (mono/di)-nucleotide at some positions. If anybody has some insight on how to do it the most efficient way(memory-wise and time-wise), then that'd be great. I use Python for my work.

May be missing something, but is there a reason you can't map the reference index to an alignment index using e.g. a list with the reference positions in the alignment? e.g. ref_posns = [3, 4, 6] if only the 4th, 5th, 7th positions are ungapped in reference; you can use this to access only those aligned positions. There is probably a more standard way to do it than that, though I don't know it off top of head. It seems rather complex to introduce an indexing scheme that doesn't have an index for many positions of the alignment. May help to know what "standard softwares" you refer to.


text:
Proteomics data Vs Transcriptomics data?
I want to use either of Proteomics or Transcriptomics data for integrating it into my kinetic model. Before proceeding, I want to know what are the advantages of using either of them so that I could make an informed decision on it!
Many studies have shown that the best we can do is integrating both transcriptomics and proteomics data with our kinetic model, but I've some time constraints and have to proceed with only one of those.
My effort and findings: I've found from discussions with researchers that gathering transcriptomics data has an amplification step which increases the chance of finding a particular one whereas gathering proteomics data has no such step but has fragmentation and then rejoining which creates many problems(such as splice variants etc) and thus leads to a loss of data. But a PostDoc told me that even after the loss of that data, I'll get more information from Proteomics data.
I want to know such type of points and want to know if these are valid or not!

You want to compare two different techniques but do you have data from both of them for your model? Theoretically they can be whatever, but then in your experiment you can get unusable data. Also it seems like you have two implementations of the model, so it might be influence the accuracy of the model. Did you try the model with simulated data of both types? (or with freely available information on the web)

I am yet to gather the data(either transcriptomics or proteomics) for my model. I have the list of genes ready for which I want it. I have only one implementation of the model only. I haven't tried the simulation with any data as I don't have it yet. I want to make an informed decision on it. Please suggest @llrs

I did some research and have put down my findings. Please add any suggestions that you want. @llrs

Ultimately, proteomics data can actually measure the actual concentrations of the enzymes that you are interested, while transcriptomics only measures their rates of change due to ongoing transcription from mRNA.
Rates of change are often treated as proxy for concentrations. This makes sense when you have a constant degradation/dilution environment and are near steady state, since then the rate of change is proportional to the converged concentration. If you're interested in more dynamic environments or believe that modulation of degradation or sequestration might play a significant role in a system, however, that assumption breaks done.
Thus, if you have to choose between proteomics and transcriptomics data, I would suggest proteomics.
As a final note, however, since the link suggests you're interested in metabolites, however, have you consider going directly to metabolomics data?  That would let you really get direct information about the particular species of interests.

This is what I found from doing some research. Comments are welcome at any point!

Capture percentage for data gathering:

Transcriptomics Data: There's an amplification step in Transcriptomics data gathering methods. Hence, it's possible to capture almost the totality of the Transcriptome using those methods(scRNA seq methods, Nanopore tech, Spatial transcriptomics, etc)
Proteomics Data: There's no amplification step in Proteomics data gathering methods. Those methods have fragmentation and a rejoining step which leads to a loss of the data.

Scalability of methods:

Transcriptomics Data: Scalable methods available
Proteomics Data: Less scalable methods for protein studies

Reference availability:

Transcriptomics Data: Fully annotated references available on consortiums such as ensembl biomart, etc
Proteomics Data: Universal and Comprehensive Human Proteome reference is still in question and because of the incomplete and inaccurate references, much of the newly generated data is being rendered useless

Uniformity

Transcriptomics Data: Lots of publications on pipelines, and their benchmarking available
Proteomics Data: Lack of uniformity across labs/research groups and lack of related literature lead to differences in protein fragmentation and solubilization and differences in algorithms to run analyses

Technical Bias

Transcriptomics Data: Many methods available to tackle cell-bias, noise in data, and batch effects, so as to capture maximum biological variability
Proteomics Data: Cell-wide Mass spectrometry has a bias towards the identification of peptides with higher concentration or contamination from other experiments

On the plus side, proteomics is actually measuring what's most relevant to your model.

@DevonRyan yes, sometimes I think that the field is focused too much on transcriptomics... I hope to see more usage (and learn more) on proteomics. On which side do you lean Akhil? BTW Did you take into account the experimental and bioinformatic experience as well as the cost of each to obtain (and store/process) the data you need? (just throwing ideas that are worth considering from a practical point of view)

I'm more inclined to use the proteomics data only but my PI has told me to go with whatever data will be of higher quality as then we can draw some inferences more confidently. Also, I didn't get your last point. Can you please elaborate on that? @llrs

Is anyone on your work (lab?) that has done RNA extractions for RNAseq or anyone that have done the preparation for mass spectrometry? Assuming you have them. How much would it cost to have the same quality of data from both techniques? Once you have the data, is there anyone that has the knowledge and resources to process it? In RNAseq you need to align to the reference (computationally expensive and time consuming), normalize/prepare it for your model while store all the data and the original files you got the fastq.gz. How do you plan to do that?

My lab is purely a computational one. We don't have the resources or people to perform the said experiments. So, my PI has asked me to gather the data from the existing literature. In the same way, we gathered the metabolites and flux data from the existing literature on Mass spectrometry and have used it to construct a kinetic model. Just to tell you I'm very well familiarised with the reference alignment, normalization, etc as I've done an internship in the scRNA seq analysis, but I don't think it'll be useful to me here in anyway. @llrs


text:
Statistical approach to compare the SNP genotyping data among set of individuals
So, I have the genotyping data of about 650,000 SNPs for 96 individuals. I already know the Y DNA haplogroup of these individuals, so to some extent, I have a gross understanding of their ancestry.
What would be the best way to go about doing this? From what I have been reading, maybe a PCA would be a good idea? Suppose individuals 1-10 belong to a particular population (unique ancestary) then their SNP genotypes should be significantly similar to each other than say individuals 11-20 who belong to another population group. So, in theory I should get two clusters in my PCA plot.
I have the genotype of each SNP for each individual. How would I go about creating a PCA plot or any other alternative way of analysis for this kind of questions.

Way better. Now one of your question is, "how to make a PCA plot". There are the plenty of tutorials out there, have you tried one? Yes? What what was the problem with it?

What is the question?

PCA sounds like a good start to see a rough population structure among your 96 individuals, this step is half a sane thing to do, but also it can give you a hint how many distinct populations there are in your dataset. I have not personally run PCA on SNP data, but it seems that there are quite lot of tutorials online. For instance this one looks  nice.
Once you know how many categories you can expect in your population, you can use Structure to classify individuals to the categories.

Sorry if I wasn't able to make it clear, but ya that's precisely what I want to do. I want to see if the individuals cluster in any manner whatsoever. Did you mean PCA by any chance in your answer? So exactly how would I go about doing it. I know what PCA does in theory but I am not able to understand how I can create a PCA plot with the data I have in hand.

@user2887 Oh, yes. I meant PCA, fixed that already. Note that you can also edit your question, to make it clearer or to add more details.

Hey, I edited the answer (added a link to a blogpost where they show how to do PCA on SNP data). Seems that quite lot of them out there if you google "pca genotype data".


text:
Question: UCSC hub: autoscale option for multiwig track
I am making a hub on UCSC Genome Browser and I am aiming to display multiple bigWig files as a single track (multiwig container). I decided to use the autoscale option because I wanted a dynamic scale: y axis changing depending on the maximum values in the current observed window.
However, I am confused as to how to the relationship between:

autoscale on and viewLimits
autoscale on and type bigWig <#> <#>

For the "type" option, the official documentation mentions that: "The remote data files of type bigWig must declare the expected signal range for the data"
So my question is, are the following assumptions correct:

If I have autoscale on, I don't need to specify viewLimits
I still have to specify "type bigWig" but you no longer have to put a range (i.e. the <#> <#>)

If this is correct, my trackDb.txt file will look as such:
track data
container multiWig
configurable on
shortLabel label data
longLabel long label data
visibility full
priority 1
type bigWig
autoScale on
aggregate transparentOverlay
windowingFunction mean
showSubtrackColorOnUi on

    track sample1
    bigDataUrl path/to/sample1.bw
    shortLabel s_sample1
    longLabel s_sample1
    type bigWig
    visibility full
    alwaysZero on
    priority 120
    color 220,180,180
    parent data on

    track sample2
    bigDataUrl path/to/sample2.bw
    shortLabel s_sample2
    longLabel s_sample2
    type bigWig
    visibility full
    alwaysZero on
    priority 120
    color 220,180,180
    parent data on

    etc

Apologies, just updated my question, I hope it's clearer.


text:
Calculating average mito.percentage for each cluster (seurat)
I have a tricky data set with cells that will have a higher percent of mitochondria genes than "typical" data sets. I would like to look at the mito percentage in each cluster without any filtering (e.g., subset(seurat_object, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mito < 10))
FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)

Done!
I would also like to calculate the average percent.mito for each cluster and place the value on the FeaturePlot/TSNE plots. I have no idea how to manipulate a Seurat object to calculate these values but I can probably figure out how to place those values on the graph -- probably with geom_text.
My Hail Mary approach is below:
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = p3, FUN = function(x) mean(p3))
CombinePlots(plots = p3)

My attempts yield the error below; I don't know if my approach is close to correct and therefore, I have spent little time working through the error message.
In mean.default(p3) : argument is not numeric or logical: returning NA
Execution halted

and
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = p3, FUN = function(x) mean(seurat_object$meta.data$seurat_clusters))
CombinePlots(plots = p3)

Error in `[[.Seurat`(x, i, drop = TRUE) :
  Cannot find 'meta.data' in this Seurat object
Calls: lapply ... lapply -> FUN -> mean -> $ -> $.Seurat -> [[ -> [[.Seurat

and
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = seurat_object, FUN = function(x) mean(seurat_object$meta.data$seurat_clusters))
CombinePlots(plots = p3)

Error in as.list.default(X) :
  no method for coercing this S4 class to a vector
Calls: lapply -> lapply -> as.list -> as.list.default```

If anyone can provide some guidance -- you will be inducted into my personal hall of fame. :)
Thank you for your time.
P.S. I am using Seurat 3.1
Update:
md = surate_object@meta.data
md_mean <- ddply(md, .(seurat_clusters), summarize, mean_mito_per=mean(percent.mito))

This was not the approach I was hoping for but I have the values -- which when expanded to the length of the original df, I can insert it into a free slot in the seurat_object.
If anyone has a better approach I would love to hear it.
UPDATE:
Using @PPK's approach
# First the cluster annotation and the tsne embeddings are merged
label.df <- cbind(as.data.frame(seurat_object@meta.data), as.data.frame(seurat_object@reductions$tsne@cell.embeddings))

# using dplyr across to calculate the mean mitochondrial percentage and
# the median tsne values per cluster
label.df <- label.df %>% dplyr::group_by(seurat_clusters) %>% dplyr::summarise(dplyr::across(percent.mito, ~ mean(.), .names = "mean_{.col}"), dplyr::across(contains("tSNE"), ~ median(.)))

# running FeaturePlot and adding the label layer
FeaturePlot(seurat_object, reduction = "tsne", features = "percent.mito") +
  geom_text(data = label.df, aes(x = tSNE_1, y = tSNE_2, label = mean_percent.mito))

Error:
* object '.col' is a function.
ℹ Input `..1` is `dplyr::across(percent.mito, ~mean(.), .names = "mean_{.col}")`.
ℹ The error occured in group 1: seurat_clusters = "0".
Backtrace:
     █
  1. └─global::initial_clustering(BM0H, "BM0H")
  2.   └─`%>%`(...)
  3.     ├─base::withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
  4.     └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
  5.       └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
  6.         └─`_fseq`(`_lhs`)
  7.           └─magrittr::freduce(value, `_function_list`)
  8.             ├─base::withVisible(function_list[[k]](value))
  9.             └─function_list[[k]](value)
 10.               ├─dplyr::summarise(...)
 11.               └

You can extract the necessary values and add them directly the plot as a second layer using plot + geom_text(). This is very similar to the inner workings of the DimPlot function with label = TRUE but allows you to use anything as label.
# I use dplyr v1.0.2 for all data frame manipulations
# First the cluster annotation and the tsne embeddings are merged
label.df <- cbind(as.data.frame(seurat_object@meta.data), as.data.frame(seurat_object@reductions$tsne@cell.embeddings))

# using dplyr across to calculate the mean mitochondrial percentage and
# the median tsne values per cluster
label.df <- label.df %>% dplyr::group_by(seurat_clusters) %>% dplyr::summarise(dplyr::across(percent.mito, ~ mean(.), .names = "mean_{.col}"), dplyr::across(contains("tSNE"), ~ median(.)))

# running FeaturePlot and adding the label layer
FeaturePlot(seurat_object, reduction = "tsne", features = "percent.mito") +
  geom_text(data = label.df, aes(x = tSNE_1, y = tSNE_2, label = mean_percent.mito))

The labels are positioned at the median position of the cells of a cluster. Typing LabelClusters (without brackets) in the R console will show you the labeling function used by Seurat that contains some additional ways of positioning the labels.
If you want to avoid overlapping labels you can use geom_text_repel from the ggrepel package.

I think this is exactly what I need @PPK. As the code is written, I get the following error: Error: unexpected ',' in: label.df <- label.df %>% group_by(seurat_clusters) %>% summarise(across(percent.mito, ~ mean(.), .names = "mean_{.col}"))," Execution halted. I believe summarise should hold both of the across items??? I have tinkered with it at bit, but have yet to resolve the error.

label.df <- label.df %>% group_by(seurat_clusters) %>% summarise(across(percent.mito, ~ mean(.), .names = "mean_{.col}"), across(contains("    tSNE"), ~ median(.)))  Yields this error:
Error: across() must only be used inside dplyr verbs.  ...                     └─dplyr::across(percent.mito, ~mean(.), .names = "mean_{.col}")` -- any tips?

Sorry for the typo in the answer. As you correctly determined the summarise function was closed too early. I modified the answer accordingly.

The second problem seems to be with some other package you are using overwriting some dplyr function. This happens pretty frequently because the functions have such generic names. I updated the answer with explicit calls to dplyr::group:by, dplyr::summarise and dplyr::across to circumvent this kind of problem. Hope this helps

There is still some issue. I will post the error at the end of my post to preserve the format.

Is your code expecting seurat_clusters to be integers rather than characters?

The code should work for both integers and characters.

I looked for solutions and it may be that the error is caused by the expressions ~ mean(.) and ~ median(.). You could try to replace this with simply mean and median.  For me both work equally well but after reading more documentation the second version may be more correct. If this solves your problem I will update the code in the answer again

Let us continue this discussion in chat.


text:
How to correctly parallelise RSeQC scripts with GNU parallel?
I have a .bam, as ouput from STAR aligner,  from which I need to extract some info using RSeQC while using all the computational resources available to increase speed. To accomplish that I am using GNU parallel. Currently I am doing it like this:

Split the main .bam per chromosome, generating files with the suffix Aligned.sortedByCoord.out.REF
bamtools split -in Aligned.sortedByCoord.out.bam  -reference

Calculate the clipping profile
ls Aligned.sortedByCoord.out.REF_* | parallel --joblog clipping_profile.log --eta --bar --max-procs 20 --jobs 20 'clipping_profile.py -i {} -s "SE" -o RSeQC.clipping.profile/clipping_profile.{.}'

This will generate a .r, .pdfand a .xlxs for each chromosome, thus ending up with a lot of files. I am forced to split the .bam per chromosome as I am not able to find a --recstart nor a --recend in the compressed .bam.
Question
Is there a better way to split a .bam file using GNU parallel and RSeQC' programs (e.g. split_bam.py), such that one gets the outputs from a single input file?
Inspired by ole.tange's thread, I have also tried this, though unsuccessfully:
cat star_outputAligned.sortedByCoord.out.bam | parallel --block 1G -k --joblog splitbam3.log --max-procs 20 --jobs 20 --pipe 'split_bam.py -i {} -r ../../My_genomes/GCA.GRCh38/GRCh38_wgEncodeGencodeAttrsV28_all_rRNA.bed12 -o splitbam3'


text:
In calculating the retention index, why do we use the character state with the lowest frequency?
To calculate the retention index for a phylogenetic tree, we use the following formula:
$$\frac{\text{maximum number of steps on tree - number of steps on the tree}}
{\text{maximum number of steps on tree - minimum number of steps in the data}}$$
To calculate the maximum number of steps, we look at the number of observed character states for a character and take the lowest number. For each character, we take the lowest number of observed traits and then sum across all characters.
Take the following dataset (columns are characters, taxa are rows):
          123456
Outgroup  000000
Taxa A    100001
Taxa B    110001
Taxa C    111001
Taxa D    111101
Taxa E    111111
Taxa F    111110

The maximum number of steps for the first character is 1, since the character state 0 is attested once and character state 1 occurs six times.
What I don't understand is why we do not use 6 for the maximum number of changes for the first character.

Could you explain what is the matrix? Each column is a feature, but what does 0 or 1 represent, a step? Perhaps I would understand better if you could show the retention index of this dataset

The retention index is used to gain insight into how well the data conform to the most parsimonious tree. 0 and 1 represent the character state.

Agreed RI is used for a parsimony tree. The minimum of steps requires searching the total tree space of the data and is not trivial to obtain without an algorithm. The maximum number of steps would computationally derrived.

found this very useful. I was wondering you could help @Kohl, how do I go about working out max steps if I have a matrix with not only 0s and 1s, but other integers such as: eg taxonA 11000203044.
taxonB 11001102034 Thanks for your help, hope you see this after it being closed for so long

The maximum number of steps (or changes) is the number of taxa with state 1 or 0, whichever is smaller. Conceptually, this is the number of character changes if each taxon evolved its state independently of the other taxa.
If for a given character the character state 0 occurred twice and character state 1 occurred five times, 2 would be the maximum number of changes.

Thank you for this. Could you just edit your answer to explain why in the scenario you sketch above the five-time occurrence of character state 1 is not used in calculating the maximum number of steps? What I find confusing is that if character state 1 occurs five times, could that not represent five independent (i.e., homoplasious) character state changes?

"maximum number of steps" corresponds to the maximum number of steps that can be reconstructed using Fitch parsimony.  Five changes from 0→1 is possible, but less parsimonious than two changes from 1→0.  Less parsimonious still, but still compatible with the data, is a million unobserved steps from 0→1 and back again.


text:
Retrieving NCBI Taxa IDs from refseq or GenBank assembly accession
I have about 10,000 genome files all named by either refseq or genbank accession number, do you know if it's possible to convert these numbers to the corresponding NCBI taxon ID or species?
for example:
GCA_000005845.2 to 79781

In the case of E.coli
Edit:
The file names look like this:
GCF_000337855.1_ASM33785v1_protein.faa.gz
GCF_000091125.1_ASM9112v1_protein.faa.gz
GCF_000184535.1_ASM18453v1_protein.faa.gz

My operating system is Ubuntu 16.04

Hi and welcome to the site. Please [edit] your question and show us a few examples of the actual file names. If you have both refseq and genbank accessions, show examples of both. Also, please mention your operating system since some solutions might depend on it.

Thanks for the edit. So, are your filenames always accession.versionassemblyprotein.faa.gz?

Yes, although they are all in different folders, for example: /refseq/archaea/GCF_000337855.1/GCF_000337855.1_ASM33785v1_protein.faa.gz

Turns out I'd already written some code that did it, human memory is a funny thing.
This uses biopython to split the field description to where the species is. May not work for all NCBI files, but seems to work on most.
import Bio
from Bio import SeqIO
from Bio import AlignIO

for record in SeqIO.parse (FILE, "fasta"):
    Speciesname = record.description.split('[', 1)[1].split(']', 1)[0]

As the question mentions python, here is a python example using biopython.
import json
from Bio import Entrez

def accession2taxid(acc: str, db="nucleotide") -> str:
    handle = Entrez.esearch(db=db, term=acc)
    record = Entrez.read(handle)
    gi = record["IdList"][0]
    handle = Entrez.esummary(db=db, id=gi, retmode="json")
    result = json.load(handle)["result"]
    taxid = result[gi]["taxid"]
    return str(taxid)

So if you want to use an assembly accession, as in your example
>>> accession2taxid("GCA_000005845.2", db="assembly")
'511145'

In your example, you said that GCA_000005845.2 should return 79781, but that Taxon ID is for "uncultured archaeon AOC672". The taxid in the example I just showed returns "Escherichia coli str. K-12 substr. MG1655", which is what you want.
Your question also mentions GenBank accessions, so here is another example
>>> accession2taxid("CU458896.1", db="nucleotide")
'561007'

You can use eutils to get this information as shown below:
esearch -db assembly -q 'GCF_000337855.1' | esummary | xtract -pattern DocumentSummary -element AssemblyAccession,Taxid
GCF_000337855.1 469382

Bumped into this question and tried @michael-hall 's answer, which works, but takes a very long time to execute if you only do one accession number per query, and there is a limit to how many you can do passing the numbers directly to esummary, so here is a version that takes a list of accession numbers, uploads them in bulk with epost, then fetches all the records at once, returning a dict mapping accession number to taxid. (I did occasionally run into random errors with the epost command, but adding some retry logic will fix it, if needed)
from Bio import Entrez

def accession_to_taxid_bulk(acc: list[str], db='nuccore') -> dict[str, str]:
    handle = Entrez.epost(db=db, id=','.join(refs))
    record = Entrez.read(handle)

    handle = Entrez.esummary(db=db,
                             query_key = record['QueryKey'], WebEnv=record['WebEnv'])
    result = Entrez.read(handle)

    return {i['AccessionVersion']: int(i['TaxId']) for i in result}


text:
pysam pileup: what reads appear in the pileup?
Cross-posted on Biostars (with no answers currently). I hope that's OK.

I thought that iff a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case.
This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.
What's wrong in my understanding?

Updated code (and updated output) moved to Gist.

Edit: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.
E00489:44:HNNVYCCXX:1:2120:26524:19135  163     5       345248  0       48M102S =       345446  350     TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG  AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)<)<--)<)--<-FAJFJFFFF-<  XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M   MD:Z:48 RG:Z:HNNVYCCXX.1       NM:i:0   MQ:i:60 AS:i:48 XS:i:48
E00489:44:HNNVYCCXX:1:2120:26524:19135  83      5       345446  60      6M2D144M        =       345248  -350    TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG  JJAF7-FJJJFJFF-AFJF7--<A--<F<7-7--FA--F-J--A-7--7-FJFJ7JFJJF<JA77-7----7----<F7-<JJF-<<FAJJJJA-7-<JAJJ<7FFJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MC:Z:48M102S    MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1        NM:i:7  MQ:i:0  AS:i:119        XS:i:69

I guess this explains why the read appears in more than 150 pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in 198 = (48 + 6 + 144) pileup columns, at each position where one of them maps?

Update: I re-ran this script, using @Bioathlete's suggestion of using the flags to distinguish the read and its mate. In case pileup() was filtering reads with low base quality, I also calculated how many positions the read has baseQ >= 13. It's 150 positions, but this read still appears in pileup columns for only 121 positions.

Can you include the raw sam lines for the read so we can determine if the issue is in pysam or your code.  That would give us truth.

@Bioathlete please see above

The two entries in the sam file represent mate pairs for the given ID.  You can tell the difference based on the sam flags.  Picard has a nice tool to determine the meaning of the flags.
You are missing coverage within the first read.  I think that it is related to the low quality bases.  samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.

Thanks for the tip about distinguishing a read and its mate by QNAME and flag—that's helpful. The read with flag 83, though, has baseQ > 13 at 150 positions, yet appears in the pileup columns for only 121 positions.

Ah but BAQ ("Base Alignment Quality"), not just regular base quality, also has to be above 13. Comparing samtools mpileup and samtools mpileup -Q0 shows that that's why reads are being filtered out.


text:
Script to allow gene set enrichment analysis of 10x genomics data in R
I have 10x single cell RNA seq data. Which R package is best suited for analysis of the 10x data matrix. What is the script to prepare the data for downstream GSEA analysis.
I have already processed samples for single cell data (for multiple other purposes). my intention is to focus on subpopulations i.e. cell type x at different time points i.e. time 0, 1, 2 etc and compare the enriched gene pathways in this one cell type at different time points. there are three files that are used (at least in the R Seurat package) with the following file extensions: .mtx, .tsv (barcodes.tsv, genes.tsv and matrix.mtx). the R Seurat package merges these files to create a "Seurat object" which helps mainly with identifying individual diffl gene expression in each barcode cell
I found a package for R: goseq that at least has been used for bulkRNA seq but I was hoping someone has used this package specifically for single cell RNA seq data, and what script has been used to merge the 3 files to create an object for analysis. pardon my terminology, but I am a novice to R and have experience primarily through package use.

Welcome to the group.  Can you specify what format you data is in and what format you need it in for the downstream GSEA.  What tools have you tried or looked into?

Hi Jay, its not clear from your question what the biological question you are trying to address with a GSEA analysis. Do you want do do enrichment of the genes that are expresesd on average in the cells of this experiement? If this is the case, how would this differ from an analysis of bulk sequencing? Do you wish to know the pathways that differ between different scRNAseq samples? Do you wish to know the pathways that differ between clusters of cells within one sample? The biological question at hand will strongly influence the downstream analysis approach.

In addition to Ian question, you mention goseq. Are you aware of the differences between pathways and gene ontologies terms? Do you have in mind any pathway database or you don't care (there are several with different focus). Are we talking about human pathways or other organisms? For other organisms there are fewer databases of pathways or less consensus. Keep updating the question. I think it can become a great question for future readers

Seruat will give you a list of genes which it thinks are upregulated in a particular cluster. Look at the functions that talk about marker genes - these functions basically do a DE analysis of the genes in one cluster compared to the others.
Then take that list and feed it to any standard GO analysis tool. Have a look at the topGO topKEGG and geneSetTest functions in the limma package, the GOStats package and the gsea. All should be suitable.  The GOSeq package is designed to compensate for gene length bias in RNA-seq experiments. As 10X only samples 3'-tags, there shouldn't be any gene-length bias in the data, so this shouldn't be an issue.

There is no purpose-built R package to perform gene set enrichment analysis on single-cell data but there does not need to be. You should be able to tools developed for bulk-RNA-Seq or microarray data, although you may not get as significant results from a sparse scRNA-Seq matrix as single-cell technologies have poor sensitivity and miss genes. What you need for a gene set enrichment analysis is:

a database of pathways (or gene sets) to test enrichment for. Packages such as GO.db or reactome.db from Bioconductor support these in R.
a gene expression matrix which you can get from the output of cellranger (as matrix.mtx). This should be normalised and corrected for batch effects. The Seurat R package is recommended to do this. Calling object@data will return the normalised seurat data matrix to pass it to other tools.
an annotation to identify the genes (depending on the reference data you used for the cellranger "transcriptome").
to define groups for differential expression analysis. You can import the clusters generated by cellranger in outs/analysis or export clusters of your choosing from the Loupe Browser as CSV files that can be imported into R. These should match the barcodes.tsv of columns of your data matrix to define the cluster for each cell barcode.
to run a differential expression analysis to compare clusters. You can use clusters to define cell groups and perform analysis with tools you would use for bulk-RNA-Seq such as limma, edgeR, or deseq. Don't forget to adjust for multiple tests. You also need to correct for batch effects between replicates and log-transform you data to ensure the results reflect biological differences. This will give you a list of differentially expressed genes.
to calculate gene set enrichment of differentially expressed genes in a gene set list such as those supported by pathway databases. You can do this like any other list of significant genes with tools used for bulk-RNA-Seq data, either with dedicated webtools or R packages such as safe or GSEA. You can also run the hypergeometric test in R yourself.

The cellranger output gives you all the data you need to use standard tools for which there are already guides on how to do differential analysis. These analyses do not differ from bulk-RNA-Seq analysis once the data has been normalised and the  clusters have been defined.
[edit:corrected typo in package name]

I realized that your question is specifically gene set enrichment analysis. However, I personally have never tried any R package doing this. Although topGO may seem to be the one that may work for you? What I have been doing is that for gene enrichment analysis using Gene Ontology, simply paste and run, or DAVID when it was still updating. Sorry, I am not an expert in incorporating biological significance study in the computational result.
It seems that you are also not quite sure what tools you should use for downstream analysis. I think the following publication may be helpful.
https://www.ncbi.nlm.nih.gov/pubmed/30228881
The same primary author also published a paper titled, Cluster Headache: Comparing Clustering Tools for 10X Single Cell Sequencing Data (For some reason I could not find the link to it, will edit the answer if I find it. )
Hope that helps.

DAVID is widely regarded as outdated. The databases that it uses are not kept up to date. This is not recommended.

@TomKelly That is absolutely right. That's why I said it was good when it was updated. Not recommended until it starts to update again.


text:
tabix indexing a GTEx variant file- ':' delimited
I have a gtex variant file- the head of which looks as follows:
phenotype_id                                    variant_id
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_13550_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14671_G_C_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14677_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_16841_G_T_b38  ...

Effectively I would like to be able to lookup the variants inside a particular gene (the file is sorted) and put them in a temporary file:
E.g. If the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt. However zgrep ENSG00000148481 file.gz is just as fast as the above...
Hence I thought tabix would be the right tool for this-
I would like to tabix index it to make lookups faster. It is gzip compressed and I will firstly have to do:
zcat gtex.txt.gz | bgzip > gtex.txt.bgz

However I am not quite sure how to proceed from there given that the data is not tab-delimited.
As a trial I tried the first 1000 lines:
zcat gtex.txt.gz | head -n 1000 | bgzip > gtex_1000.txt.bgz

./tabix -p bed gtex_1000.gz #index as a bed file
[get_intv] the following line cannot be parsed and skipped: chr1:15947:16607:clu_36198:ENSG00000227232.5        chr1_13550_G_A_b38     ......

./tabix -p vcf gtex_1000.gz #index as a vcf file
Indexing as a bed file results in a warning while indexing as a vcf gives no warning yet either way when I try to retrieve a sequence:

./tabix test.gz chr1:15000:17000

It returns nothing.
I am starting to think that I will just have to write a script that splits on the ':' and writes the data to a new file.... and then index that file- which will take a huge amount of time
... Does anyone know of a trick to index the files with unconventional delimiting?

Suppose that you could successfully creat an index file for your txt file, which program would be used to read it?

@PhoenixMu see the way that I have used sed in the updated question... this was my plan... although this was the same speed as grep

I think tabix only works for bed and vcf format. Your file is not in bed format?

Yes it's not in bed format as the coordinate columns are : delimited and the file is too big to justify producing a copy with properly delimited columns- I am in the process of making an index file myself and will post an answer if it works

I proceeded as follows:
import gzip
import sys

qtls=['sqtls','eqtls']
tissues=['Adipose_Subcutaneous', 'Adipose_Visceral_Omentum', 'Adrenal_Gland', 'Artery_Aorta', 'Artery_Coronary', 'Artery_Tibial', 'Brain_Amygdala', 'Brain_Anterior_cingulate_cortex_BA24', 'Brain_Caudate_basal_ganglia', 'Brain_Cerebellar_Hemisphere', 'Brain_Cerebellum', 'Brain_Cortex', 'Brain_Frontal_Cortex_BA9', 'Brain_Hippocampus', 'Brain_Hypothalamus', 'Brain_Nucleus_accumbens_basal_ganglia', 'Brain_Putamen_basal_ganglia', 'Brain_Spinal_cord_cervical_c-1', 'Brain_Substantia_nigra', 'Breast_Mammary_Tissue', 'Cells_Cultured_fibroblasts', 'Cells_EBV-transformed_lymphocytes', 'Colon_Sigmoid', 'Colon_Transverse', 'Esophagus_Gastroesophageal_Junction', 'Esophagus_Mucosa', 'Esophagus_Muscularis', 'Heart_Atrial_Appendage', 'Heart_Left_Ventricle', 'Kidney_Cortex', 'Liver', 'Lung', 'Minor_Salivary_Gland', 'Muscle_Skeletal', 'Nerve_Tibial', 'Ovary', 'Pancreas', 'Pituitary', 'Prostate', 'Skin_Not_Sun_Exposed_Suprapubic', 'Skin_Sun_Exposed_Lower_leg', 'Small_Intestine_Terminal_Ileum', 'Spleen', 'Stomach', 'Testis', 'Thyroid', 'Uterus', 'Vagina', 'Whole_Blood']

gene_old=''

def index(tissue,qtl,gene_old):
        counter_old,counter=1,1
        if qtl=='sqtls':
                path='/path/to/file/'+tissue+'.v8.sqtl_allpairs.txt.gz'
        else:
                path='/path/to/file/'+tissue+'.allpairs.txt.gz'
        f = gzip.open(path, 'r')
        next(f)
        for line in f:
                try:
                        line=line.decode()
                        line=line[:line.find('\t')]
                        gene=line[line.rfind('ENS'):].split('.')[0]
                        #print(gene,gene_old,counter_old,counter)
                        with open(qtl+'_'+tissue+'.idx','a') as indexed:
                                 indexed.write(str(gene_old)+'\t'+str(counter_old)+'\t'+str(counter)+'\n')
                        gene_old=gene
                        counter +=1
                        counter_old=counter
                except IndexError:
                        print('tissue,qtl,gene_old,counter')
                        counter+=1
        f.close()

for qtl in qtls:
        for tissue in tissues:
               qtl=row.qtls
               tissue=row.tissues
               print(qtl, tissue)
               index(tissue,qtl,gene_old)

The above script writes an index for each gene which should look as follows (took 2 days to finish all files):
ensid    row_start    row_end
......
ENSG00000230337     1526391 1534329
ENSG00000171819 1534330 1542321
ENSG00000198793    1542322 1550178
ENSG00000120942  1550179 1558010
.....

And to get a gene of interest you will have to look up the coordinates for an ensemblid/gene and for example:
zcat Brain_Nucleus_accumbens_basal_ganglia.v8.sqtl_allpairs.txt.gz | sed -n '1550179,1558010p;1558011q' > ENSG00000120942.txt

Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?

Hi @Colin thanks for getting to me- I just want to be able to take out all variants inside a particular gene- for example if the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt
Hence I thought tabix would be the right one for this

Here is an example of a tabix file that has been hacked up to allow "free text query". http://target.wustl.edu/ENCFF000ARO_ALL.iteres.loci.gz  So they can query tabix file.txt "A(n)" and because there are many lines starting with A(n) it will return all of them

Can also check out other options for text indexing https://unix.stackexchange.com/questions/3086/command-line-friendly-full-text-indexing


text:
Canu assembly not making a single consensus?
I've downloaded reads from this BioProject. Using canu with default parameters (no correction), I've got 4 contigs, none of which really look like the reference plasmid here.
The command I used was:
canu -p ip40a -d ip40a_assembly -useGrid=false -maxThreads=6 -genomeSize=175k -nanopore-raw reads.fastq

Granted, the original publication for this sequence used minimap/miniasm but my questions are:
A. Why do I get multiple contigs? Doesn't canu generate a single assembled sequence?
B. This is a plasmid but canu says suggestedCircular=false. What's going on here?

Looking that the link you provided only 30% of the reads appear to be e. coli.  I would suggest filtering the reads for those that align to e. coli and see where that gets you.  Did the original publication describe the method they used for data cleaning?
A:
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR3385327

Unidentified reads: 68.82%
Identified reads: 31.18%
cellular organisms: 30%
Bacteria: 30%
Proteobacteria: 28.09%
Gammaproteobacteria: 19.7%
Enterobacterales: 13.3%
Enterobacteriaceae: 12.54%
Vibrionales: 0.01%
Viruses: 1.17%

B:
Without a good assembly canu won't be able to determine if the assembled genome is circular.  I would ignore this until you have a good example.

Sadly the original publication only seems to say "assembled by miniasm/minimap" (along those lines) and I'm trying to figure out the pipeline as we speak -- there's no indication on whether the reads were error corrected, etc.

Also, how would I filter the E. coli-identified reads?

I would align the reads to the plasmid reference and see how many align and what the coverage of the reference looks like.  Are there any regions of the reference that are uncovered.  how even is the coverage across the reference.  These are things that can help you understand what to expect from the assembly.

It seems like it's actually the opposite that @thestatnoob needs to do: exclude E. coli reads, instead of selecting for them

Correct he has the plasmid reference he is trying to assemble.  He needs to select for those reads assuming the alignment coverage shows an assembly is even worthwhile.

The E. coli genome is closer to 4.6 Mb in size; if there are a large quantity of host reads in there, then they'll take over the reads and there won't be enough for your target plasmid. By default, it will only try to assemble from a maximum of ~40X coverage of the specified genome size (higher values do not add any more information for assembly).
It's possible that miniasm will take all reads into account, rather than just a portion of reads. You may need to filter out the E. coli reads for this assembly to work properly.
For what it's worth, Canu's default parameters do include a correction step.

Yep, but the reads aren't for an E. coli genome. It's for a plasmid within E. coli, which is definitely within the 170-180k region (in fact, I got the number from the ENA...). I'm basically just trying to reconstruct the consensus plasmid that the authors have, using canu instead of miniasm.

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


text:
How to trim adapter sequences from GSE65360 in order to map the reads?
I am trying to map single-cell chromatin accessibility data from doi:10.1038/nature14590, obtained using scATAC-seq, to the reference genome. Example paired-end reads are here.
The reads contain various adapter sequences, which the authors "trimmed from FASTQs using custom python scripts" [ref]. I am having problems reproducing this trimming. They used Nextera Transposase Adapters

TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
  GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG

and custom barcode adapters which are listed [here]. An excerpt from this list:

Custom Barcodes Adapter 1 (index i5):
  AATGATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT
  AATGATACGGCGACCACCGAGATCTACACTAACCAAGTCGTCGGCAGCGTCAGATGTGTAT
  AATGATACGGCGACCACCGAGATCTACACAGGTTGCCTCGTCGGCAGCGTCAGATGTGTAT
  Custom Barcodes Adapter 2 (index i7):
  CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG
  CAAGCAGAAGACGGCATACGAGATCCCAACCTGTCTCGTGGGCTCGGAGATGTG
  CAAGCAGAAGACGGCATACGAGATCACCACACGTCTCGTGGGCTCGGAGATGTG

Analyzing the FASTA-files with FastQC showed that the reverse complement of those sequences was over-represented in the data and also confirmed Nextera Transposase Adapters.
I tried to trim the sequences using cutadapt:
cutadapt --minimum-length=20 -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATA -A ATACACATCTGACGCTGCCGACGA -A CTGTCTCTTATA -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz SRR1780164_1.fastq.gz SRR1780164_2.fastq.gz
with

CACATCTCCGAGCCCACGAGAC being the reverse complement of the 2nd fix part of the i7 adapters
ATACACATCTGACGCTGCCGACGA being the reverse complement of the 2nd fix part of the i5 adapters
CTGTCTCTTATA for the Nextera Transposase Adapters

FastQC does not indicate overrepresented sequences in the resulting FASTA-files, however still indicates Nextera Transposase Adapters in the _1 files. On the other hand, after TrimGalore removes CTGTCTCTTATA, Nextra Adapters are not indicated anymore, for some reason.
Furthermore, mapping the reads using bowtie2
bowtie2 -x GRCh38_primary_assembly_index -1 SRR1780164_1_trimmed.fastq.gz -2 SRR1780164_2_trimmed.fastq.gz > SRR1780164.sam
results in a low alignment rate
1352 (2.79%) aligned concordantly exactly 1 time
3420 (7.07%) aligned concordantly >1 times
11.25% overall alignment rate

which gives me the impression I did something wrong during adapter trimming.
Any help would be much appreciated.
--
Edit: As I am not sure about the sequences to use for trimming, I also tried -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA, as found in a suggestion for Nextera Adapters, however it did not remove the sequence over-representation in the -1 file. -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA solves adapter-issues according to FastQC, but also does not give good mapping results.

First the Nextera adapters and the custom barcode adapters overlap each other
Nextera 1                                     TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Barcode 1   GATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT

Nextera 2                                    GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Barcode 2    CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG

Looking at some read 1s and read 2s you were correct that the adapters appear at the 3' end of the read reverse complimented (The |> points to the adapter in the sequence).  You can see that CTGTCTCTTATACACA would be the ideal sequence to use for the 3' end of the reads.  Since the Nextera adapter and the barcode overlap you only need to use the one that is adjacent to the insert.
SRR1780164_2.fastq: GCCCGAGGGGACGCTGCCCCACGTGCTCCAGG |> CTGTCTCTTATACACATCTGACGCTGCCGACGATGCACGAAGTG
SRR1780164_1.fastq: GGGCTGGACAGCGGCGGGGACGCGGACCTGCAG |> CTGTCTCTTATACACATCTCCGAGCCCACGAGACACAAAGTGA

Of the 64k reads there are ~11k with the Nextera adapter at the 3' end based on a rough sequence grep
On the 5' end of the read there is some Nextera adapter present but it is not reverse complimented and there are only a handful of reads with the adapter present ~250 total between read 1 and read 2.
SRR1780164_2.fastq: TGGGCTCGGAGATGTGTATAAGAGACAG <| CTCCCCCGCCTCCAGCATCCGGGCGAGGTAGTGCATCGACGCGTCCAC

I used this command:
cutadapt --minimum-length=20 -A CTGTCTCTTATACACA -a CTGTCTCTTATACACA -G GATGTGTATAAGAGACAG
-g GATGTGTATAAGAGACAG -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz
SRR1780164_1.fastq SRR1780164_2.fastq

Which gives a trimming result of
   === Summary ===

Total read pairs processed:             64,832
  Read 1 with adapter:                  16,721 (25.8%)
  Read 2 with adapter:                  16,996 (26.2%)
Pairs that were too short:                  46 (0.1%)
Pairs written (passing filters):        64,786 (99.9%)

I think this best represents the library construct and the trimming of the adapter sequences

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


text:
Is there any tools that can generating FBA model from genome sequences or proteome sequences?
Suppose I have all sequences of genes or proteins of certain species, are there any free tools I can use to generate a flux balance analysis model out of them?
What I was thinking is more like to annotate all the enzymes based on some annotation method. Then mapping these enzymes to known meta-metabolic networks that we know for instance the meta-metabolic network in KEGG database.

To generate flux balance, one should now the stoichiometric numbers of all the reactions that are taking place, as well as the conditions (that could affect the reactions). I am not sure we have enough knowledge, time and computationally power to do this.
At the moment I'm only aware of a paper simulating the Mycoplasma genitalium cell. But I'm not sure it was worth the effort as I haven't seen replicated or redone with other species, perhaps the model was not accurate enough to be worth it.
The tools used are here. The licenses are not a standard, but at least they seem free.

Thanks. What I was thinking is more like to annotate all the enzymes based on some annotation method. Then mapping these enzymes to known meta-metabolic networks that we know for instance the meta-metabolic network in KEGG database.

Which would be your input? A list of genes? This is trivial to do, but you don't know how many active sites are in an enzyme nor if there is any inhibitor, and you would need to know the initial concentrations to see if you reach a steady state or not, and what is your end point is also a variable you should define


text:
10X Illumina demultiplexing sample sheet issue
Also posted on biostars.
I am trying to use cellranger or bcl2fastq to convert the .bcl files that I got from single cell analysis run into fastq files for further analysis. I needed to generate sample_sheet.csv and so I used the following tool:
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/bcl2fastq-direct
which I copy-pasted into a text file using vim. sample_sheet.csv looks like that:
[Header]
EMFileVersion,4

[Reads]
26
8
98

[Data]
Lane,Sample_ID,Sample_Name,index,Sample_Project
1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406
3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180406
4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180406
4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180406
4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180406
4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180406
4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180406
4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180406
4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180406

Then I ran cellranger using bash script:
#!/bin/bash

 FLOWCELL_DIR="/scratch/nv4e/kipnis/180403_NB501830_0158_AHN3LLBGX5"
 OUTPUT_DIR="/scratch/nv4e/kipnis/fastq"
 SAMPLE_SHEET_PATH="/scratch/nv4e/kipnis/sample_sheet_2.csv"

cellranger mkfastq --id="AHN3LLBGX5" --run=${FLOWCELL_DIR} --csv=${SAMPLE_SHEET_PATH} --output-dir=${OUTPUT_DIR}

It fails with the following output:

_stderr file shows:

and that my sample_sheet.csv can not be parsed:

It could be because CLRF needs to be converted into LF, but I tried two commands to do this unsuccessfully:
tr -d '\r' < input > output and perl -pi -e 's/\r\n/\n/g' input
From the following thread:
https://unix.stackexchange.com/questions/277217/how-to-install-dos2unix-on-linux-without-root-access
The above commands successfully strip the \r from the line endings, but the error from bcl2fastq remains.
What am I doing wrong? Any suggestions would be greatly appreciated.

Update

I looked into the folder that is generated when the cellranger runs and found two .csv files generated by the cellranger. One of them, samplesheet.csv, contains these ^M characters that are shown in the error:

Update

Maybe one of the issues was that I was submitting regular .sh file on the SLURM cluster. However, I wrote the right .slurm file and submitted it yesterday with sbatch command:

However, today I got the message in my email:

And neither _stderr, nor generate_fastq_id.err or generate_fastq_id.out files contain anything. I thought initially that I set not enough time and it timed out, however, it is not true: I set the time to 24h and it has been running for just 13h before failing.

Update

I searched over all of the files in the generated folder with the find linux command and actually found the error files that show the same error as before: issue with the sample_sheet.csv file. This error file is AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u179acba3d8/_stderr.
So, right now the issue seems to be the following: I am feeding in the right sample_sheet.csv without the ^M symbols, but cellranger transforms it into one that contains ^M and then uses it giving the error. Above I see also which: no configureBclToFastq.pl and ERROR: bcl2fastq::common::Exception: 2018-Apr-09 13:33:14: Inappropriate ioctl for device, so maybe my input parameters are just wrong to cellranger and it has nothing to do with the sample_sheet.csv?

Update

I tried to use --sample-sheet instead of --csv since I am feeding in complex sample sheet. No luck.

Update

Here you can find all the error document in text:

which: no configureBclToFastq.pl in (/sfs/lustre/scratch/nv4e/cellranger/samtools_new/1.6:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-fastq/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-ref/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c9/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/bin:/sfs/lustre/scratch/nv4e/cellranger/lz4/v1.8.0:/sfs/lustre/scratch/nv4e/cellranger/martian-cs/2.3.2/bin:/sfs/lustre/scratch/nv4e/cellranger/STAR/5dda596:/sfs/lustre/scratch/nv4e/cellranger/samtools_new/1.6:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-fastq/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/cellranger-tiny-ref/1.2.0:/sfs/lustre/scratch/nv4e/cellranger/miniconda-cr-cs/4.3.21-miniconda-cr-cs-c9/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/lib/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/tenkit/bin:/sfs/lustre/scratch/nv4e/cellranger/cellranger-cs/2.1.1/bin:/sfs/lustre/scratch/nv4e/cellranger/lz4/v1.8.0:/sfs/lustre/scratch/nv4e/cellranger/martian-cs/2.3.2/bin:/sfs/lustre/scratch/nv4e/cellranger/STAR/5dda596:/scratch/nv4e/cellranger:/scratch/nv4e/bcl2fastq/build/bin:/scratch/nv4e/spark/bin:/scratch/nv4e/scala/bin:/sfs/nfs/blue/nv4e/private/bin:/sfs/nfs/blue/nv4e/anaconda2/bin:/sfs/nfs/blue/nv4e/.local/bin:/usr/lib64/qt-3.3/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/opt/slurm/current/bin:/opt/slurm/current/sbin:/opt/singularity/current/bin:/opt/rci/bin:/opt/rci/sbin:/opt/nhc/current/sbin:/share/rci_apps/common/bin:/share/resources/HPCtools/)
  BCL to FASTQ file converter
  bcl2fastq v2.20.0.422
  Copyright (c) 2007-2017 Illumina, Inc.
  2018-04-09 18:31:29 [7fcdb39f37c0] Command-line invocation: bcl2fastq --minimum-trimmed-read-length 8 --mask-short-adapter-reads 8 --create-fastq-for-index-reads --ignore-missing-positions --ignore-missing-filter --ignore-missing-bcls --use-bases-mask=Y26,I8,Y98 -R /scratch/nv4e/kipnis/180403_NB501830_0158_AHN3LLBGX5 --output-dir=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u21f1cbe9bf/files/fastq_path --interop-dir=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/BCL2FASTQ_WITH_SAMPLESHEET/fork0/chnk0-u21f1cbe9bf/files/interop_path --sample-sheet=/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/PREPARE_SAMPLESHEET/fork0/chnk0-u09e7cbe891/files/samplesheet.csv -p 6 -r 6 -w 6
  2018-04-09 18:31:29 [7fcdb39f37c0] INFO: Minimum log level: INFO
  2018-04-09 18:31:29 [7fcdb39f37c0] INFO: Sample sheet: '/scratch/nv4e/kipnis/AHN3LLBGX5/MAKE_FASTQS_CS/MAKE_FASTQS/PREPARE_SAMPLESHEET/fork0/chnk0-u09e7cbe891/files/samplesheet.csv'
  2018-04-09 18:31:29 [7fcdb39f37c0] ERROR: bcl2fastq::common::Exception: 2018-Apr-09 18:31:29: Inappropriate ioctl for device (25): /scratch/nv4e/bcl2fastq/src/cxx/include/common/CsvGrammar.hpp(92): Throw in function bcl2fastq::common::CsvGrammarAttribute bcl2fastq::common::parseCsvData(Iterator, Iterator) [with Iterator = __gnu_cxx::__normal_iterator >; bcl2fastq::common::CsvGrammarAttribute = std::vector > >]
  Dynamic exception type: boost::exception_detail::clone_impl
  std::exception::what: Could not parse the CSV stream text:
   ^M
  ^M
  [Reads]^M
  26^M
  98^M
   ^M
  ^M
  [Data]^M
  Lane,Sample_ID,Sample_Name,index,Sample_Project,Original_Sample_ID^M
  1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M
  2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M
  3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M
  4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409,SI-GA-B4_1^M
  4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409,SI-GA-B4_2^M
  4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409,SI-GA-B4_3^M
  4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409,SI-GA-B4_4^M
  4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409,SI-GA-B5_1^M
  4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409,SI-GA-B5_2^M
  4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409,SI-GA-B5_3^M
  4,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409,SI-GA-B5_4^M

What happens if you remove the duplicate emtpy lines (just keep one before each section)?

There are no duplicate lines in the input file, but inside the generated _stderr file they are somehow present

Is that the output from _stderr after removing \r? Because the ^M characters indicate a line break problem. If they go away after running the tr/perl command, then there's probably an additional problem

Yes, it is after fixing the sample sheet.

Did you manage to install the configureBclToFastq.pl in the path as reported in the _stderr? BTW, if you make your error message as text it would be easier for us to search and quote specific parts of the error: CTRL+SHIFT+C and CTRL+SHIFT+V in terminal

No, I have not yet tried installing it because I do not know how. Could you give any links of how to do that and where to get the file?

There's no reason to install the old bcl2fastq (the configureBclToFastq.pl that's mentioned). Can you directly run bcl2fastq with your sample sheet (i.e., without using cellranger)?

So, we figured it out. If you delete everything until the [Data] tag in the sample sheet, leaving the [Data] tag, it starts working. Here is the full sample_sheet.csv that is working:
[Data]
Lane,Sample_ID,Sample_Name,index,Sample_Project
1,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
1,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
1,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
1,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
1,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
1,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
1,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
1,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
2,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
2,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
2,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
2,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
2,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
2,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
2,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
2,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
3,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
3,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
3,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
3,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
3,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
3,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
3,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
3,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409
4,SI-GA-B4_1,17R,ACTTCATA,Chromium_20180409
4,SI-GA-B4_2,17R,GAGATGAC,Chromium_20180409
4,SI-GA-B4_3,17R,TGCCGTGG,Chromium_20180409
4,SI-GA-B4_4,17R,CTAGACCT,Chromium_20180409
4,SI-GA-B5_1,19RL,AATAATGG,Chromium_20180409
4,SI-GA-B5_2,19RL,CCAGGGCA,Chromium_20180409
4,SI-GA-B5_3,19RL,TGCCTCAT,Chromium_20180409
4,SI-GA-B5_4,19RL,GTGTCATC,Chromium_20180409

It is very unfortunate and disappointing that 10X website is generating wrong sample sheet. Still we got very low mapping rates and trying to figure out what went wrong.

Update

I was supplying --sample-sheet= property, but it should be --samplesheet=. Now it is running fine (except for very low mapping rates) with all the headers on.

Update

It turned out that we were given the wrong indexes, so we made a sample sheet using sample sheet generator again and even with --samplesheet= correct parameter it was not running. Convertion of CLRF with perl or tr did not help. However, concise format of the sample sheet worked out well which looks like the following:
Lane,Sample,Index
1-4,17R,SI-GA-B8
1-4,19RL,SI-GA-B9

It needs to be supplied with --csv= parameter.

I found that if you copy and paste the 10X sheet directly from their website it inserts some odd characters on some of the blank lines.  I usually edit the file in vim and delete them.

Related to @GWW, I suspect it was the  ^M lines (the ones with a space and then line end). Either way, I suggest that you email 10x support and report this issue. It should be easy for them to fix.

It looks like it's a weird unicode character "U+00A0" / NO-BREAK SPACE.

That's weird. I copied and pasted from their website and never had issues. But I also run mac2unix/dos2unix on every sample sheet, so maybe that is a sufficient fix.

Apart from supporting only Windows-style and Linux-style line endings and a restricted character set as mentioned above, the sample sheet has to be in the correct encoding which is UTF-8 (or its subset US-ASCII).
Many command lines support these lookup and conversion commands:
# lookup
file -i my_sample_sheet.csv

# conversion
iconv -f ISO-8859-15 -t UTF-8 in.csv > out.csv


text:
tx2gene file for tximport issue
I am trying to use tximport to get raw count matrices of a bulk RNA-seq sequencing experiment. I followed the tutorial:
https://bioconductor.org/packages/3.7/bioc/vignettes/tximport/inst/doc/tximport.html
Where they generate tx2gene in the following way:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

So, I reasoned to do it the same way. But since my data is mouse one, I used a different gene annotation file:
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

The transcript annotation file that I used when running salmon is gencode.vM13.transcripts.fa.sidx. Then when I try getting txi object:
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

It raises the error:

Seems like I am creating tx2gene in the wrong way. Maybe I need to use a different library, not TxDb.Mmusculus.UCSC.mm10.knownGene. Any suggestions would be greatly appreciated.

Update

I looked into the quant.sf files that I got from salmon and it looked the following way:

A colleague adviced me to strip away everything that follows the dot in the Name column, so it started to look:

However, still no luck. The identifiers in the k <- keys(txdb, keytype = "TXNAME") look like that:
"uc007afg.1" "uc007afh.1" "uc007afi.2" "uc011wht.1" "uc011whu.1"
"uc057aty.1"

I suppose I need to convert these identifiers first to ENS-like ones, but do not know how to.

The solution was to use a different library as I expected. It should be library(TxDb.Hsapiens.UCSC.hg19.ensGene). And also, as specified in the question, to remove everything after the dot in the Name column of the quant.sf files.

Do your ids in the expression file (from gencode) match up with the gene IDs in the transcript file (from UCSC?) I'm willing to bet that they don't. Stay consistent throughout your process. Either pull down the gencode transcript to gene translations, or do your pseudoalignments to UCSC transcripts.

I already got the salmon aligned data, that was aligned to the gencode file. Now I do not know what TxDb library file I need to use.

see answer here


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:
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.