Bacterial genome annotation of a clinical isolate strain? So I'm basically so new to this that I'm just trying to find out what tools, methods, and keywords I should go look up by myself.

I have a unique strain of a bacteria.

I was given RNAseq data for this unique strain.

I want to analyze the RNAseq data, but I need an annotated genome to do that.

As a resource, I do have the annotated genome of the common lab strain.

The common strain should be pretty similar to the unique strain, but I don't know how similar.

Can anyone tell me what tools/methods might be applicable so I can go look them up?

For a quick (but reliable) analysis, I'd recommend using Kallisto or Salmon to quantify isoform read counts using the transcriptome of the lab strain.

If you have a concern about transcripts that are in your sample but not in the lab strain, you can do two Trinity assemblies: one fully de-novo, and one genome-guided. These transcripts can then be used as a reference assembly for Kallisto/Salmon, and counts compared to the lab-strain results. Note that the Trinity assemblies will only apply to the particular environments that are in your samples; if a transcript isn't expressed, it won't be in the assembly.

Any transcripts in the de-novo assembly that are not in the genome-guided assembly are potentially novel to your strain. However, care must be taken in interpreting this: if there is any other contamination in the sample, those contaminant transcripts will also be assembled.

I'm not too familiar with annotation pipelines, but based on a couple of Twitter posts it looks like Prokka might be a reasonable start. Torsten Seeman (lead developer of Prokka) has done a post about alternative annotation pipelines:

http://thegenomefactory.blogspot.co.nz/2013/03/bacterial-genome-annotation-systems.html


Phylogenetics software that can represent collapsed nodes as triangles I am trying to represent a large (100-300+ sequences) phylogenetic tree in a way that gets my point across. I figured that reducing the number of branches (by collapsing clades which include highly similar sequences) would be a good idea. The software that I'm currently using (iTOL) technically has an option for that ('collapse node'), but I'm not exactly a fan of the particular graphical representation, so I'd like to have a look at other options. Thank you for your time!

For Mac Figtree here is highly recommended, it is also available for Windows and Linux here. The technique is called "collapsing clades" and is the only package I know that does it. I would recommend using it, because I would assume it is averaging the branch length of a given clade and the 'base' of the triangle is proportional to the number of taxa. Thus it recognised as being very reliable, albeit the precise algorhithm it uses ... I don't know :-O

With FigTree you get two different "triangle options" in the menu, which give different shaped triangles for collapsing the number of lineages within the select clade you chose. The big one, is the algorhithm appropriately described above. It has many options for producing publication quality diagrams, albeit I then import the output into a graphics program for further editing (there are shortcomings IMO in colouring and fonts. Then again for free is better than good.

If you are outputing a tree within a script, ETE3 (Python) will do a "collapsing of clade" option, but it isn't pretty triangle stuff. ETE3 will produce a publication quality tree diagram from this and reduce to 'low resolution', which is very useful in large analyses of loads of trees.

In summary, I don't think FigTree has any competition in the context you are seeking.

Moreover, you can configure your bin, or alias file so you can automatically open a treefile just typing,

figtree treefile.tre

Also you can reduce you font size to next point 1 or 2 (okay you can do that in anything), if you are attempting optical analysis on loads of taxa. Its also really stable. I just find it really useful, like a "jazz standard" its a "phylogeny standard".


A GTF annotation for GRCH375d version I have been given a bunch of .bam files, that come from the alignment on GRCh37d5 of human genome by Star tool. I want to extract raw read counts from these bam files by FEATURECOUNTS tool so I need a GTF. I googled a lot but I did not find GTF for GRCh37d5. I used another GTFs and I obtained weird features. Please one of you help me in finding this specific GTF. I converted a bam to fastq and I used Salmon to quantify gene expression by a transcripome.fasta but rate of mapping was 0.06

The most recent annotation of GRCh37 assembly by NCBI RefSeq is here: https://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/GRCh37.p13_interim_annotation/

You can download the annotation in GFF3 format from that path and use it with featurecounts tool by first converting it to the SAF format using the commands shown below:

$ zgrep -v '^#' annotation.gff.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($3=="exon"){print $1,$4,$5,$7,$9}' \
    | sed -r 's/ID.*GeneID:([0-9]*).*/\1/g' \
    | awk 'BEGIN{FS="\t";OFS="\t"}{print $5,$1,$2,$3,$4}'
    > annotation.saf

You should confirm that the chromosome sequence identifiers in the BAM file and the SAF file match. NCBI annotation has identifiers in the accession.version format whereas UCSC uses chr# format and Ensembl uses #. For this, you can use the -A parameter from featureCounts. Assuming that your BAM file uses UCSC-style identifiers, you can generate a 'chromosome name alias file' that can be used with featureCounts as follows:

## get NCBI assembly report for GRCh37
wget -nv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.24_GRCh37.p12/GCF_000001405.24_GRCh37.p12_assembly_report.txt

## generate chr_alias.csv file
grep -v '^#' GCF_000001405.24_GRCh37.p12_assembly_report.txt \
    | dos2unix \
    | awk 'BEGIN{FS="\t";OFS=","}($10!="na"){print $10,$7}' \
    > chr_alias.csv


How to decide number of neighbors and resolution for Louvain clustering I am using Louvain clustering (1,2) to cluster cells in scRNAseq data, as implemented by scanpy.

One of the parameter required for this kind of clustering is the number of neighbors used to construct the neighborhood graph of cells (docs).

Larger values result in a more global view of the manifold, leading to lower number of clusters, while reducing the number of neighbors goes in the opposite direction. However, it is unclear how to choose this parameter.

The resolution parameter seems to work in the opposite way.

Do you know of any methodology and/or rule-of-thumb to define these parameters? E.g. depending on the size of the dataset?

Levine, Jacob H., et al. "Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis." Cell 162.1 (2015): 184-197.
Blondel, Vincent D., et al. "Fast unfolding of communities in large networks." Journal of statistical mechanics: theory and experiment 2008.10 (2008): P10008.

A general rule of thumb is that in order to improve the variance $n$ times you need $n^2$ neighbours. This is only applicable if you consider the $n^2$ nearest neighbours of a cell to be biologically identical (i.e. "similar enough"); if your data includes 10 types of cells with 10 cells each, then using the 20 nearest neighbours for smoothing will obscure the data.

As far as I know, there is no single best answer to this question. I would suggest trying different numbers and sticking to what agrees more with the biology of the dataset.


Set new Idents based on gene expression in Seurat and mix n match identities to compare using FindAllMarkers I am relatively new to Bioinformatics and scRNA-seq data analysis. I am using Seurat V3 to analyze a scRNA-seq dataset in R. Currently, I have merged three scRNA-seq samples from the same donor into one Seurat object, All_Samples. We'll call them:

Uninfected
Virus1
Virus2

All_Samples <- merge(x = Uninfected, y = c(Virus1, Virus2))

Uninfected cells did not receive virus (i.e. express no viral genes), the other two samples were infected with a virus introducing new genes into the cells. I would like to accomplish two things.

First: I want to create a new Ident for Virus1 and Virus2 samples based on expression of a viral gene  of interest: "GeneA". All cells from Virus1 and Virus2 samples that have "GeneA" expression > 0.5 would be labeled "Pos", those with "GeneA" < 0.5 would be labeled "Neg" in a  the new Column called GeneA.

I am able to subset the objects based on GeneA expression that applies to all samples in the object. For example:

All_Samples_GeneA_Pos <- subset(All_Samples, subset = GeneA > 0.5)
All_Samples_GeneA_Neg <- subset(All_Samples, subset = GeneA < 0.5)

I assume I would have to modify All_Samples@meta.data based on this post and this post, but I admit I'm not entirely sure how to implement the suggested answers into my data.

Second: Once I have separated the data from Virus1 and Virus2 into "Pos" and "Neg" cells for "GeneA", I want to look for differentially expressed genes between all Uninfected cells and "Pos" cells from Virus1 or Virus2. I have been using something like this to compare samples:

 FindMarkers(object = All_Samples, group.by = 'virus', ident.1 = "1", ident.2 = "2")

But how would I do something equivalent to this:

 FindMarkers(object = All_Samples, ident.1 = "Uninfected", ident.2 = "Virus1_GeneA_Pos")

Or maybe

 FindMarkers(object = All_Samples, ident.1 = "Uninfected_GeneA_Neg", ident.2 = "Virus1_GeneA_Pos")

Let me know if I can clarify any points. I appreciate any help!

For the first question, you can use ifelse() to create a new column in the meta.data slot:

All_Samples@meta.data$new_column <-
  ifelse(
    rownames(All_Samples@meta.data) %in% colnames(All_Samples_GeneA_Pos),
    "GeneA_Pos",
    "GeneB_Pos"
  )

colnames(seurat_object) provides a vector of cell names in a given Seurat object. Here whatever cell that is in the All_Samples_GeneA_Pos object would be GeneA_Pos and whatever is not GeneB_Pos. To better control the behavior, you can use a "nested" ifelse(); you can have another ifelse() instead of the "GeneB_Pos" bit above.

For the second question, the rationale is the same as above, you can define an extra column with the "identities of interest", then use these identities with SetIdent() and then in turn using two identities of interest in the differential expression call.


Do any publicly available databases detail protein structure and functional domains? I am interested in finding a database that takes a gene or protein name as input (possibly with the option to specify transcript) and gives information about the protein's functional domains in terms of either specific residue ranges, genomic coordinates, etc. Ideally, the database would include 3D conformational information as well as the functional domains. Does such a database exist? It seems that this information can be inferred for a given protein using various publications and references but I have yet to find a single collection containing these data all in one place.

Some of this information (at least some domains, active sites, etc) is available from UniProt.

If you want to download their whole database, you can search without specifying any terms and then click the Download button.


Retrieve detailed gene descriptions Given a list of gene IDs, how do you retrieve the gene description, summary and other detailed information in R?

Description of how to use biomaRt here. Filter by your list of gene IDs. Get the descriptions etc as attributes. Eg:

> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> IDs <- c("BRCA2","BRAF")
genedesc <- getBM(attributes=c('external_gene_name','description'), filters = 'external_gene_name', values = IDs, mart =ensembl)
> genedesc
  external_gene_name                                                                      description
1               BRAF B-Raf proto-oncogene, serine/threonine kinase [Source:HGNC Symbol;Acc:HGNC:1097]
2              BRCA2                  BRCA2, DNA repair associated [Source:HGNC Symbol;Acc:HGNC:1101]


Differentiating molecules based on peptide sequence? How to annotate? I want to differentiate between classical class I and non classical class I MHC molecules in a model organism using well conserved structural features within classical MHC I molecules (eg intradomain disulfide bridges within class I a1 domain). To do this I need to align the peptide sequences of transcripts of that model organism with transcripts of classical MHC I human genes and see which genes from those transcripts have those structural features. So I am looking for a tool to align and annotate those structural features in a replicable manner
Additionally if you have any other bioinformatic ideas /computational analyses on how to find classical MHC Class I molecules on an organism which is poorly annotated that would highly appreciated!

Knowing almost nothing about the biology, I think that the most straightforward way to do this is:

train profile HMMs on your different classes of protein sequences (in human maybe, but ideally including MHCs from a large number of organisms)
translate your transcripts into proteins and run them against your trained profile HMMs
confirm that the profiles adequately discriminate the two classes and adjust if needed.

For an example of how to do this (on MHC), you can see here.
These approaches are then easily extensible to further organisms, as you have the profile HMMs.
It's possible that a simpler approach would consist of just finding hits against human MHC and then just checking for the existence of sequence patterns corresponding to those structural motifs, but I doubt that would work well in any significantly divergent organism.


SNP-phenotype association analysis only using the SNPs on a specific gene I know that the GWAS association p-value threshold is 1e-8. This makes sense because in GWAS, you make ~1 million hypothesis (i.e. use that many SNPs in the association test). However, let's say I do not do genome-wide analysis but only use the 100 SNPs from a specific gene to check for association with my phenotype. Then, do you think a raw p-value of 1e-4 (which becomes 1e-2 after multiple hypothesis correction) is a reliable p-value to come to a conclusion that there is an association between that SNP and the checked phenotype? Obviously, if I checked all 1 million SNPs instead, this p-value would be insignificant after correction, but I am not checking 1 million SNPs, as I mentioned.

I know that the GWAS association p-value threshold is 1e-8

This may be a common threshold of statistical significance that is used, but it's definitely not an absolute value. It's a hack to try to work around many issues with GWAS associated with testing millions of SNPs. Unfortunately, the most relevant issue in GWAS (for spurious significance) is unexpected shared genetic structure in the cases or controls, and this cannot be excluded by a p-value threshold, regardless of how low it is set.

In your 100-SNP example, you are correct that the bonferroni correction would be a threshold of 1e-4 (for a desired significance threshold of 0.01), but it's important to be cautious about GWAS results, even when the p-value indicates otherwise. Trust no one, especially yourself. P-values should never be used as an indication of importance (which suggests that most GWAS results should probably be reconsidered), and are better used in conjunction with other evidence for identifying relevant / significant tests.

My recommendation for GWAS (or any similar bulk genotyping study, as in your 100-SNP subset) is that at least a bootstrap sub-sampling process is carried out, to make sure that observations found within the comparison of interest are at least well replicated when comparing smaller sub-groups of the cases and controls. More information on that can be found in the poster I presented at Queenstown Research Week last year, and detail about the implementation in my draft preprint (which will probably never be properly published, given how old the research is).


lower mapping rates in salmon v0.13 compared to previous versions Hi there :) Thanks for the tool!
I recently updated to the new salmon (from 0.8... its been a couple years) and I noticed that my mapping percentages change dramatically between the two versions.
For example, using the default settings in v0.8, I see a mapping rate of 96.25% however upon upgrading to v0.13, these mapping rates go down to 82-84% depending on which of the parameters I play with (I've permutated  the validateMappings, incompatPrior, consensusSlack, minScoreFraction, ma, etc). Below I have a snippet of one of these tests.  I'm currently assembling the unmapped reads to see if I can determine why these reads were not mapped.
What do you think is the most likely cause for this dramatic change in performance?
Perhaps it could be the transcript index? (I used a kmer of 31 and my reads are 100 bp).
UPDATE:
I have run the indexing and quasi-mapping  using salmon v0.8, v0.10 (when validateMappings was released), v0.12, v0.13 and had mapping rates of 96, 90, 90 and 83 respectively.  So it looks like a gradual decay of mapping rates.  Or alternatively, a gradual reduction in spurious mapping?
Seems like what ever improvements from 0.12 to 0.13 changed the mapping rates significantly.  I will go ahead with both mappings (v0.12 and v0.13) to see if there are differences in the DESEQ2 outputs.
Thanks a lot!
(while I have presented values here for one set of fastq files, the same pattern is seen throughout all 72 of my samples from different experiments)
for F in $FILES ; do
  R1=${F}_L001_R1_001.fastq.gz
  R2=${F}_L001_R2_001.fastq.gz
  salmon quant \
      -p 35 \
      -i $ref_dir/$transcript_index \
      -l A \
      -1 $R1 \
      -2 $R2 \
      -o $(basename $F).VM.rFB4.IP9.quant \
      --incompatPrior 9.9999999999999995e-21 \  # default of v0.8
      --writeUnmappedNames \
      --validateMappings \
      --rangeFactorizationBins 4

done

Update salmon is now released with a --allowDovetail option in  0.13.1 and inspected via a -z dump.

[Most likely answer based on comments attached to the question]
Hi Courtney, We've looked into the data and identified the source of the different mapping rates. Specifically, the cause is discordant reads; in this case reads that are dovetailing with respect to their positions on the reference. If you run Bowtie2 with "RNA-seq" flags (i.e. disallowing discordant mappings), then you get essentially the same mapping rate as salmon here. Is allowing dovetailing reads something that should be very important in your analysis? If so, we can consider what would be necessary to add back support, but generally, they provide evidence of tenuous quality.


Base abbreviations other than ACTGU in sequence file Short question:
where can I find what F and J and < and - and 7 mean in a Sequence file?
Longer version:
I was looking at this sequence file I got from NCBI (it's an sra file? if that is the correct terminology) and I used the sratoolkit to print out a few spots
looking at them I expected to see ACGU (this was meant to be an RNA sequence), but what I see are like

AAFFFJJJFJJJJJJJJJJJJJJJJJJJJJFJJF<JJ77AJJJJ-JJJJJFF7FA--F7FAJJ7<AJ-FJA<J<JJFA<<77F-JJJJJJF<JF<7-FF-A-----A-77JJ7FJJF-FFFAJ<<-FJJFAJAFFAFFFFF<7A7JFJJ7AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJFJJJJJJJJJJJJAA<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJAAFFJJFJJJJJJJJJJJJJJFJ-FJF

Where can I find a reference to help me "understand" this?
I mean I assume - means missed or unknown? or is it meant to be already aligned and these are deletions?
I appreciate any pointer to get started. as you must have observed I am a total noob.
I got the file from: https://www.ncbi.nlm.nih.gov/sra/SRX8156420[accn]
Then I did a
 ./fastq-dump -X 5 -Z ERR4557575

Follow-up question:

Where do you look for the quality scores ranges and such, I guess as @winnie2k suggested I should convert this to FASTQ and make my life easier, but still wanted to understand the file better.

I usually convert SRA files to FASTQ before doing any further processing, but the tooling for converting SRA files is cumbersome.
To simplify things, I recommend you download the FASTQ file directly. The ENA currently has a nicer interface for doing these things. Here is a link to download the FASTQ files for the SRA sequencing run associated with the SRA experiment you are describing: https://www.ebi.ac.uk/ena/browser/view/SRR11588968
A description of the FASTQ format can be found here: https://en.wikipedia.org/wiki/FASTQ_format


probeset to probeset mappings between Affymetrix arrays I am interested in identifying mappings between different types of Affymetrix arrays. I am aware that mappings between gene and probeset can be extracted using Ensembl's Biomart database.

Ensembe gene id ENSG00000181019 maps to
1. AFFY HG-U133_Plus_2's 210519_s_at
2. AFFY HuGene-1_0-st-v1's 8002303

Is there any way to extract mappings of probeset ids between two arrays(eg: HG-U133_Plus_2, HuGene-1_0-st-v1)?

eg: 210519_s_at and 8002303

If your question is: can probeset IDs from different platforms be mapped to one another in a similar way as mapping probesets to genes, then the answer is: Yes. BioMart allows you to map almost anything that has an ID to anything else that has an ID.

You can use BioMart either via the web interface or programatically. A brief guide to using the web interface, for mapping HG U133 Plus 2 to HuGene 1.0:

Go to the start page
Select H. sapiens Ensembl genes for your database; Human genes (GRCh38.p10) for your dataset
Click Filters in the left-hand column
Expand "REGION", scroll down to GENE and select "Input microarray probes/probesets ID list [Max 500 advised]"
Select AFFY HG U133 PLUS 2 probe ID(s) and either copy/paste or upload a list, one per line
Click Attributes in the left-hand column
Scroll through, uncheck what you don't want to see and choose what you do, for example Gene Stable ID, AFFY HuGene 1 0 st v1 probe and AFFY HG U133 Plus 2 probe
Finally click "Results" in the menu at the top of the left-hand column

And you should see a result like this (you'll need to click the "Results" button).

You should not expect that there be a one-to-one mapping, for two reasons:

Genes have multiple transcripts and probesets are mapped to each transcript
Some transcripts have more than one probeset: in this case, the HuGene IDs 8002301 and 8002303 map to transcripts for this gene

Finally: here's a programmatic example using R/BioMart:

library(biomaRt)
mart.hs <- useMart("ensembl", "hsapiens_gene_ensembl")
results <- getBM(attributes = c("ensembl_gene_id", "affy_hugene_1_0_st_v1", "affy_hg_u133_plus_2"),
                 filters = "affy_hg_u133_plus_2",
                 values = c("210519_s_at"),
                 mart = mart.hs)

results
  ensembl_gene_id affy_hugene_1_0_st_v1 affy_hg_u133_plus_2
1 ENSG00000181019               8002301         210519_s_at
2 ENSG00000181019               8002303         210519_s_at


Genome assembly of SRR12196449 with SPAdes I am trying to assemble the run SRR12196449 with SPAdes. The description of their project is:

This project expected to standardize a method for amplification and
sequencing of FIV genome in a simple way, allowing a broader analysis
to increase the knowledge on the biology and evolution of the virus
and virus-host interaction.

This is from School of Veterinary Medicine; University of Sao Paulo.
I have tried in three ways:

Using a fasta file downloaded from NCBI (using Filtered Download method)

SPAdes-3.14.0/spades.py    --only-assembler -s sra_data.fa -o raw-fiv1

Using a fastq file also downloaded from NCBI using Filtered Download method.

SPAdes-3.14.0/spades.py   --12 sra_data.fastq.gz -o raw-fiv1-b

Downloading the original file and dumping to a fastq using sratools.

SPAdes-3.14.0/spades.py    --12 SRR12196449.fastq -o raw-fiv1-c

In the third case I get a very bad assembly, likely because it needs trimming. In the two former cases, I get an assembly with a max contig length of around 4k-5k bps. But when I use quast to evaluate the assembly I obtain a largest alignment of ~500bps. I use this as reference.
Furthermore, if make a blast search it finds that the 4k-5k contig matches other FIV sequences at ~90% . In between these other sequences, there is at least one (MF370550.1) submitted by the University of Sao Paulo (likely another run from another sample in the same project).
On the other hand, if I assemble the run using another software I obtain a contig that is an almost perfect match to my reference sequence (~99%).
I have a few questions.

Am I using SPAdes correctly? Is there a better way to use it?
Is it possible that the sequences that blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match (not because these are correct)?
What other option could explain these circumstances that I could be missing?

EDIT:
Some clarifications regarding some comments from @MaximilianPress. I can confirm that I used the same reference genome in both cases. I used that command to get the quast results:
quast-5.0.2/quast.py -r raw-fiv1/sequence.fasta SPAdes-3.14.0/raw-fiv1/contigs.fasta

The other assembler is an overlap-layout-consensus algorithm that I implemented. The methodology is similar: I use exactly the same input fasta file.
I have also tried using the SPAdes --plasmid flag and that doesn't improve much the result. Max contig length is now 3642 bps and largest alignment 737 bps, covering only 26% of the reference.
BTW, I made also tests with other virus sequences. I am just interested to know if I am using SPAdes correctly so that I can compare the performance of my algorithm Here there is a link to additional cases with other genomes. There you can find all the data and results used (raw files, assembly results and quast results)
Another clarification: I have also got perfect assemblies using SPAdes for other runs corresponding to other viruses... I am following exactly the same procedure in this case.

Update 2:
I looked into this a little more, with the various data sources.
This is related in part to the answer submitted by OP juanjo75es, in addition to discussion on chat. I don't entirely understand the logic, but the general thrust seems to be that SPAdes makes weird assemblies somehow.
Some notes that I made:
REFERENCE ASSEMBLIES

FIV sequence U11820.1 was deposited in 1996, before SPAdes existed. Unclear assembly method.
FIV sequence MN630242 was deposited in 2020, but was assembled with CLC workbench.
These two sequences align together ok with minimap2. Low exact identity, it's true. This is the PAF output:

MN630242.1      8977    194     8964    +       U11820.1        9462    575     9326    1130    8778    60      tp:A:P  cm:i:116        s1:i:1125
       s2:i:55 dv:f:0.1761

These two assemblies are syntenic, and both have similar gene calls using prokka. GFFs:

MN630242.1      Prodigal:2.6    CDS     253     1605    .       +       0       ID=AMGANPBD_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00001;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     1656    4868    .       +       0       ID=AMGANPBD_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=AMGANPBD_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
MN630242.1      Prodigal:2.6    CDS     4861    5616    .       +       0       ID=AMGANPBD_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00003;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     5891    8461    .       +       0       ID=AMGANPBD_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00004;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     8626    8790    .       +       0       ID=AMGANPBD_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00005;product=hypothetical protein

U11820.1        Prodigal:2.6    CDS     634     1983    .       +       0       ID=EANIPDKN_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00001;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     1995    5246    .       +       0       ID=EANIPDKN_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=EANIPDKN_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
U11820.1        Prodigal:2.6    CDS     5239    5994    .       +       0       ID=EANIPDKN_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00003;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     6269    8830    .       +       0       ID=EANIPDKN_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00004;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     8904    9152    .       +       0       ID=EANIPDKN_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00005;product=hypothetical protein

I can share FAA files of the proteins if needed.
READ ASSEMBLY
I also assembled the indicated reads using SPAdes. For reference it is a ~9Kbp virus, but this is a 4.3Mbp assembly. There is a lot of non-virus sequence in there. The second-largest contig is a shuttle vector. the third-largest contig is cat (host). Many more are cat, so I think it's pretty cat-oriented. The original authors used CLC workbench, so I guess that just worked a lot better in this instance, even in the presence of all the contamination. Unclear why, it appears that CLC works similarly to OP's assembly tool.
I aligned these to MN630242.1 with minimap2. If I understand OP, they are unhappy about the mapping of these contigs to this genome reference. The identities are mostly high (though there is indeed some oddness with the largest contig, which only finds very low coverage- at the same time that there are overlapping contigs with very high ID?), and they cover the entire reference genome:
MN630242.1  8977    2155    3199    -   NODE_14_length_1054_cov_3786.620280 1054    2   1046    1021    1044    60  tp:A:P  cm:i:189    s1:i:1021   s2:i:378    dv:f:0.0031
MN630242.1  8977    7893    8626    -   NODE_25_length_743_cov_3814.258117  743 8   739 713 733 60  tp:A:P  cm:i:129    s1:i:713    s2:i:318    dv:f:0.0045
MN630242.1  8977    72  4904    -   NODE_1_length_4942_cov_25.814123    4942    34  4863    628 4832    41  tp:A:P  cm:i:7s1:i:628  s2:i:512    dv:f:0.1703
MN630242.1  8977    3376    3897    +   NODE_39_length_526_cov_5254.155388  526 3   524 512 521 0   tp:A:S  cm:i:8s1:i:512  dv:f:0.0023
MN630242.1  8977    6270    6673    +   NODE_82_length_409_cov_1.744681 409 3   406 357 403 0   tp:A:P  cm:i:53 s1:i:357    s2:i:351    dv:f:0.0185
MN630242.1  8977    6852    7223    -   NODE_87_length_381_cov_23.334646    381 6   376 355 371 2   tp:A:P  cm:i:5s1:i:355  s2:i:351    dv:f:0.0079
MN630242.1  8977    6187    6558    -   NODE_88_length_380_cov_2338.675889  380 7   378 351 371 0   tp:A:S  cm:i:5s1:i:351  dv:f:0.0085
MN630242.1  8977    6852    7308    +   NODE_56_length_471_cov_0.933140 471 15  471 351 456 0   tp:A:S  cm:i:43 s1:i:351    dv:f:0.0397
MN630242.1  8977    7332    7707    +   NODE_86_length_389_cov_3506.202290  389 5   380 334 375 17  tp:A:P  cm:i:5s1:i:334  s2:i:309    dv:f:0.0111
MN630242.1  8977    5181    5568    -   NODE_60_length_464_cov_1.005935 464 50  437 314 387 17  tp:A:P  cm:i:41 s1:i:314    s2:i:282    dv:f:0.0357
MN630242.1  8977    7277    7707    +   NODE_76_length_438_cov_0.919614 438 5   435 309 430 0   tp:A:S  cm:i:40 s1:i:309    dv:f:0.0437
MN630242.1  8977    6953    7285    +   NODE_94_length_345_cov_2.903670 345 8   340 301 332 0   tp:A:S  cm:i:50 s1:i:301    dv:f:0.0110
MN630242.1  8977    5744    6056    -   NODE_113_length_316_cov_2.169312    316 4   316 289 312 35  tp:A:P  cm:i:4s1:i:289  s2:i:244    dv:f:0.0096
MN630242.1  8977    5615    5927    +   NODE_108_length_322_cov_1103.400000 322 6   318 283 312 58  tp:A:P  cm:i:4s1:i:283  s2:i:209    dv:f:0.0143
MN630242.1  8977    6543    6813    -   NODE_147_length_279_cov_2592.519737 279 8   278 261 270 21  tp:A:P  cm:i:4s1:i:261  s2:i:236    dv:f:0.0032
MN630242.1  8977    8544    8872    +   NODE_100_length_333_cov_2259.189320 333 1   329 252 328 0   tp:A:P  cm:i:3s1:i:252  s2:i:245    dv:f:0.0318
MN630242.1  8977    7601    7853    -   NODE_200_length_255_cov_25.125000   255 0   252 250 252 12  tp:A:P  cm:i:3s1:i:250  s2:i:236    dv:f:0.0017
MN630242.1  8977    6648    6968    -   NODE_106_length_324_cov_2.538071    324 4   324 248 320 0   tp:A:P  cm:i:3s1:i:248  s2:i:248    dv:f:0.0303
MN630242.1  8977    5940    6193    -   NODE_162_length_266_cov_1809.287770 266 6   259 246 253 5   tp:A:P  cm:i:4s1:i:246  s2:i:240    dv:f:0.0033
MN630242.1  8977    5416    5664    -   NODE_264_length_255_cov_9.531250    255 5   253 243 248 0   tp:A:P  cm:i:4s1:i:243  s2:i:242    dv:f:0.0033
MN630242.1  8977    7221    7466    -   NODE_247_length_255_cov_12.539062   255 1   246 241 245 0   tp:A:P  cm:i:4s1:i:241  s2:i:241    dv:f:0.0032
MN630242.1  8977    4830    5071    -   NODE_1121_length_245_cov_11.711864  245 4   245 239 241 0   tp:A:P  cm:i:4s1:i:239  s2:i:235    dv:f:0.0042
MN630242.1  8977    7769    8013    +   NODE_677_length_251_cov_3.951613    251 6   249 238 244 7   tp:A:P  cm:i:4s1:i:238  s2:i:230    dv:f:0.0075
MN630242.1  8977    4901    5148    +   NODE_270_length_255_cov_8.726562    255 8   255 238 247 0   tp:A:P  cm:i:4s1:i:238  s2:i:234    dv:f:0.0073
MN630242.1  8977    8705    8948    -   NODE_382_length_254_cov_5.763780    254 7   250 234 243 15  tp:A:P  cm:i:4s1:i:234  s2:i:217    dv:f:0.0062
MN630242.1  8977    6065    6307    +   NODE_410_length_254_cov_3.488189    254 1   243 227 242 10  tp:A:P  cm:i:3s1:i:227  s2:i:210    dv:f:0.0122
MN630242.1  8977    5062    5239    -   NODE_7003_length_180_cov_9.207547   180 3   180 172 177 17  tp:A:P  cm:i:2s1:i:172  s2:i:157    dv:f:0.0066
MN630242.1  8977    11  136 +   NODE_21321_length_128_cov_1195.000000   128 2   127 125 125 3   tp:A:P  cm:i:2s1:i:125  s2:i:123    dv:f:0

Only a (relatively) few contigs align well (figure), as might be expected with heavy contamination, and of those some are overlapping:

To look into the odd behavior of the biggest contig, I ran prokka on it too and found the expected genes for its position, especially the largest gene in the virus, dut (Deoxyuridine 5'-triphosphate nucleotidohydrolase). I then took the protein sequences of the 3 assemblies and aligned them with clustalo. They're all quite similar:
CLUSTAL O(1.2.4) multiple sequence alignment

NODE1_KJHFFCBH_00001         ------MEKRPEIQIFVNGHPIKFLLDTGADITILNRKEFIIGNSIENGKQNMIGVGGGK
MN630242_AMGANPBD_00002      -------------MIFVNGYPIKFLLDTGADITILNRRDFQVKNSIENGRQNMIGVGGGK
U11820.1_EANIPDKN_00002      MGTTTTLERRLEIQIFVNGHPIKFLLDTGADITILNRKDFQIGNSIENGKQNMIGVGGGK
                                           *****:*****************::* : ******:**********

NODE1_KJHFFCBH_00001         RGTNYINVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
MN630242_AMGANPBD_00002      RGTNYTNVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISDKI
U11820.1_EANIPDKN_00002      RGTNYINVHLEIRDENYRMQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
                             ***** ***********: **************************************:**

NODE1_KJHFFCBH_00001         PIVKVRMKDPTQGPQVKQWPLSNEKIEALTEIVERLEQEGKVKRADPNNPWNTPVFAIKK
MN630242_AMGANPBD_00002      PIVKVKMKDPNKGPQIKQWPLSNEKIEALTEIVERLEKEGKVKRADPNNPWNTPVFAIKK
U11820.1_EANIPDKN_00002      PIVKVRMRDPIQGPQVKQWPLSNEKIEALTDIVERLESEGKVKRADPNNPWNTPVFAIKK
                             *****:*:** :***:**************:******.**********************

NODE1_KJHFFCBH_00001         KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQWKKQVTVLDIGDAYFTIPLDPDYAP
MN630242_AMGANPBD_00002      KSGKWRMLIDFRELNKLTEKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
U11820.1_EANIPDKN_00002      KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
                             ************ *****:**************** ************************

NODE1_KJHFFCBH_00001         YTAFTLPRKNNAGPGKRYVWCSLPQGWVLSPLIYQSTLDNILQPYIKQNPELDIYQYMDD
MN630242_AMGANPBD_00002      YTAFTLPRKNNAGPGRRYIWCSLPQGWILSPLIYQSTLDNIIQPFIRQNPQLDIYQYMDD
U11820.1_EANIPDKN_00002      YTAFTLPRKNNAGPGRRYVWCSLPQGWVLSPLIYQSTLDNILQPFIKQNSELDIYQYMDD
                             ***************:**:********:*************:**:*:** :*********

NODE1_KJHFFCBH_00001         IYIGSNLSRQEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
MN630242_AMGANPBD_00002      IYIGSNLNKKEHKEKVGELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWTIQQKQL
U11820.1_EANIPDKN_00002      IYIGSNLSKKEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
                             *******.::***:** ************************************:******

NODE1_KJHFFCBH_00001         EIPERPTLNDLQKLAGKINWASQTIPKLSIKALTHMMRGDQKLDSIREWTEEAKKEVQKA
MN630242_AMGANPBD_00002      DIPEQPTLNELQKLAGKINWASQAIPDLSIKALTNMMRGNQNLNSIREWTKEARLEVQKA
U11820.1_EANIPDKN_00002      EIPERPTLNELQKLAGKINWASQTIPDLSIKELTNMMRGDQKLDSIREWTVEAKREVQKA
                             :***:****:*************:**.**** **:****:*:*:****** **: *****

NODE1_KJHFFCBH_00001         KEAIEKQAQLYYYDPNRELYAKISLVGPHQLCYQVYHKNPEQILWYGKMNRQKKRAENTC
MN630242_AMGANPBD_00002      KKAIEGQVQLGYYDPSKELYAKLSLVGPHQISYQVYQRNPERILWYGKMSRQKKKAENTC
U11820.1_EANIPDKN_00002      KEAIEKQAQLNYYDPNRGLYAKLSLVGPHQICYQVYQKNPEHILWYGKMNRQKKKAENTC
                             *:*** *.** ****.: ****:*******:.****::***:*******.****:*****

NODE1_KJHFFCBH_00001         DIALRACYKIREESIVRIGKEPVYEIPASREAWESNIIRSPYLKASPPEVEFIHAALNIK
MN630242_AMGANPBD_00002      DIALRACYKIREESIIRIGKEPKYEIPTSREAWESNLINSPYLKAPPPEVEYIHAALNIK
U11820.1_EANIPDKN_00002      DIALRACYKIREESIIRIGKEPMYEIPASREAWESNLIRSPYLKAPPPEVEFIHAALSIK
                             ***************:****** ****:********:*.****** *****:*****.**

NODE1_KJHFFCBH_00001         RALSMVQEVPILGAETWYIDGGRRQGKAARAAYWTNTGRWQVMEIEGSNQKAEVQALLMA
MN630242_AMGANPBD_00002      RALSMIKDAPILGAETWYIDGGRKLGKAAKAAYWTDTGKWQVMELEGSNQKAEIQALLLA
U11820.1_EANIPDKN_00002      RALSMIQDAPITGAETWYIDGSRKQGKAARAAYWTDTGKWQIMEIEGSNQKAEVQALLLA
                             *****:::.** *********.*: ****:*****:**:**:**:********:****:*

NODE1_KJHFFCBH_00001         LKEGPEEMNIITDSQYILNIMNQQPDLMEGIWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
MN630242_AMGANPBD_00002      LKAGPEEMNIITDSQYVINIILQQPDMMEGIWQEVLEELEKKTAIFIDWVPGHKGIPGNE
U11820.1_EANIPDKN_00002      LKAGSEEMNIITDSQYILNIINQQPDLMEGLWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
                             ** * ***********::**: ****:***:*******:*** *****************

NODE1_KJHFFCBH_00001         EVDKLCQTMMIIEGDGILNKRPEDAGYDLLAAQEIHILPGEVRIVPTRTRIMLPKGYWGL
MN630242_AMGANPBD_00002      EVDKLCQTMMIIEGDGILDKRSEDAGYDLLAAKEMHLLPGEVKVIPTGVKIMLPKGYWGL
U11820.1_EANIPDKN_00002      EVDKLCQTMMIIEGEGILDKRSEDAGYDLLAAQETHFLPGEVRIVPTKTRIMLPKGHWGL
                             **************:***:** **********:* *:*****:::** .:******:***

NODE1_KJHFFCBH_00001         IMGKSSIGNKGMDVLGGVIDEGYRGEIGVIMINLSKKSTTILEKQKVAQLIILPCKHESL
MN630242_AMGANPBD_00002      IIGKSSIGSKGLDVLGGVIDEGYRGEIGVIMINLSRKSITLLEQQKIAQLIILPCKHEVL
U11820.1_EANIPDKN_00002      IMGKSSIGSKGMDVLGGVIDEGYRGELGVIMINLTKKSITILEKQKVAQLIILPCRHESL
                             *:******.**:**************:*******::** *:**:**:********:** *

NODE1_KJHFFCBH_00001         EQGEIIMNSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNIPRIVAEEI
MN630242_AMGANPBD_00002      EQGKVIMDSERGDKGYGSTGVFSSWVDRIEEAEINHEKFHSDPQYLRTEFNLPKMVAEEI
U11820.1_EANIPDKN_00002      QQGEIQMDSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNLPRIVAEEI
                             :**:: *:****:**:*****************:*****************:*::*****

NODE1_KJHFFCBH_00001         KRKCPLCRIRGEQVEGKLKIGPGIWQMDCTHFNGKIIIVAIHVESGLLWAQIIPQETADC
MN630242_AMGANPBD_00002      RRKCPVCRIRGEQVGGQLKIGPGIWQMDCTHFDGKIILVAIHVESGYIWAQIISQETADC
U11820.1_EANIPDKN_00002      KRKCPLCRIRGEQVGGQLKIGPGIWQMDCTHFNGKIIIVAVHVESGFLWAQIIPQETAEC
                             :****:******** *:***************:****:**:***** :***** ****:*

NODE1_KJHFFCBH_00001         TLKAIMQLVSTHNVTEIQTDNGPNFKNQKIEGLLSYMGIKHKLGIPGNPQSQALVENANS
MN630242_AMGANPBD_00002      TVKAVLQLLSAHNVTELQTDNGPNFKNQKMEGALNYMGVKHKFGIPGNPQSQALVENVNQ
U11820.1_EANIPDKN_00002      TVKALLQLICAHNVTELQTDNGPNFKNQKMEGLLNYMGIKHKLGIPGNPQSQALVENANN
                             *:**::**:.:*****:************:** *.***:***:**************.*.

NODE1_KJHFFCBH_00001         TLKVWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYTQQESLRIQDYFSQLP
MN630242_AMGANPBD_00002      TLKAWIQKFLPETTSLENALALAVHCLNFKQRGRIGGMAPYELLAQQESLRIQEFFSKIP
U11820.1_EANIPDKN_00002      TLKAWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYIQQESLRIQDYFSQIP
                             ***.************:******:**********:* ******  ********::**::*

NODE1_KJHFFCBH_00001         QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
MN630242_AMGANPBD_00002      QKLQAQWIYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIKRVPEPCALPE
U11820.1_EANIPDKN_00002      QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
                             ***  **:*****************************************:******:***

NODE1_KJHFFCBH_00001         GDE
MN630242_AMGANPBD_00002      GDE
U11820.1_EANIPDKN_00002      GDE
                             ***

Overall, the level of homology is quite high at the protein level. There's maybe a little evidence that this contig's protein is closer to U11820.1 than to MN630242 (though there are also several positions where U11820.1 is the outgroup). So I'm not sure why pieces of this contig are so hard to align at the DNA level without further investigation, but honestly overall of these are looking like very similar viral sequences, as might be expected.
I also aligned the big contig to both U11820.1 and MN630242 at the DNA level in 3-way clustalo alignment it is ok. Not sure what the difference is there. There is no obvious reason to prefer one reference over the other, according to my eye. When I have clustalo output a clustering solution the big contig is essentially equidistant between the two references, so I don't buy that U11820.1 is "better" as a reference:
(
MN630242.1:0.116151
,
(
NODE1_rc:0.0504856
,
U11820.1:0.0504856
):0.0656651
)
;

As I noted in the comments, it's unfortunate that QUAST and SPAdes aren't working very well. I don't know what the deal is there. But it seems like overall the assembly is quite close to both references. I don't have any intuition for why SPAdes is doing what it's doing, or why alignment at the DNA level is acting weird (lots of artificial gaps?). But I don't think it's a pathology of SPAdes particularly, it seems like it's just something weird that happened with this dataset interacting with DNA aligners. Possibly something went strangely with read deposition?
Update:
See OP's self-answer as well (and also discussion in the comments).

Based on the (extremely informative!) follow-up updates that you have provided, I think that we can tentatively answer your questions:
The tl;dr is that you can check your reads to make sure they aren't weird and you can check your output sequence to see if it's weird.

Am I using SPAdes correctly? Is there a better way to use it?

I think that you are using SPAdes correctly, based on your experience getting ~finished viral genomes in the past with the same workflow. (THis is also my experience using SPAdes on viral genomes with Illumina PE data, is that it "just works".)
It is possible that you are using some other tools in your workflow in a suboptimal way. For example, maybe your trimming is not quite correct (it does seem based on your answer that you are already trimming). The question in my mind is then whether there is some non-intuitive thing in the reads that should still be trimmed.
I would suggest running FASTQC on your reads for this run and also your reads for the other runs that have worked, to see if something weird is going on. FASTQC directly checks for adapter sequence, weird GC composition, quality scores etc. That will give you a lot of information.

Is it possible that the sequences that blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match (not because these are correct)?

This is formally possible, of course. Our main way of knowing that you're getting the right genome is... does it look like other genomes that have been generated in the same way?
If you are really worried about this, I would suggest directly inspecting the sequence to see what the matches are. If you like, you can then blast the matches to see if they
Another approach would be to annotate the genome and make sure that it "looks like" a phage genome. E.g. it has the expected genes from FIV. Prokka is very easy to run and works ok on phage.
A slightly self-serving suggestion is to use PhageTerm to check your reads against your final genomes to ensure that they behave like phage- e.g. they have termini in the correct places, etc. (it does this with test coverage.) I am a contributor to PhageTerm, so I happen to know it a bit. Likely there are other options.

What other option could explain these circumstances that I could be missing?

As I've suggested, I think that that leaves the reads. How sure are you that they are comparable to your other runs? They look like they should be sufficient from my glance at SRA, but maybe someone bumped the MiSeq while it was running. I've given suggestions above for checking with FASTQC.


Access workdir defined on command line from within Snakefile Snakemake provides access to a workflow object within a Snakefile. This allows one to, for example, have dynamic programmatic access to the directory containing the Snakefile (via the workflow.basedir attribute). Is there a similar way to access the working directory that is specified on the command line using --directory?

Through a bit of exploration I was able to figure out that the workflow object has a workdir() method, but it appears this is used internally by Snakemake for registering working directories, not for access.

UPDATE: It looks like the working directory is stored internally as the _workdir attribute of the workflow object. However, this attribute is unset (default value None) when I'm creating static/global objects in the Snakefile. Is it impossible to access the working directory prior to building the workflow DAG?

You could query Python for the working directory within the Snakefile, since evidently Snakemake changes the actual working directory of the process to the one specified by --directory.    For example, using a stub Snakefile with just these two lines:

import os
print("Current working directory: " + os.getcwd())

And on the command line:

$ snakemake --directory working-directory
Creating specified working directory working-directory.
Current working directory: /home/jesse/dev/snakemake-wd/working-directory
Building DAG of jobs...
Nothing to be done.
Complete log: /home/jesse/dev/snakemake-wd/working-directory/.snakemake/log/2019-08-02T142009.349338.snakemake.log


hisat2 --rna-strandness option and downstream htseq-count analysis I've got some doubts on the hisat2 --rna-strandness option and its output for downstream analysis. Please see below.

I understand that the --rna-strandness option produces an XS tag to indicate where a transcript is from (on the + or - strand) for downstream transcriptome assembly analysis. I have a paired-end stranded sequencing library that was aligned to the genome using hisat2 without specifying the --rna-strandness (in other words, the default unstranded was the usage). Following this, the reads were assigned to genes using htseq-count and this time "-s reverse" was specified given the strand-specific sequencing assay type.

Would the above handling affect the alignment and counting results given the default usage of --rna-strandness in hisat2 followed by htseq-count -s reverse on a strand-specific assay? Since --rna-strandness is for transcriptome assembly using the XS tags generated and htseq does not use XS tags for counting, I presume there should be no practical impact from the above. Could you also shed light on this? in case I may have been overlooking other facts of the usages of the tools.

To help verify the above, I re-aligned and counted the reads from 2 samples by switching on --rna-strandness RF in hisat2. I attach the alignment and count features info. below for assessment.

Overall alignment rate of Sample 1: 94.52% (--rna-strandness RF) vs.94.12% (--rna-strandness unstranded)
Overall alignment rate of Sample 2: 94.57% (--rna-strandness RF) vs.94.15% (--rna-strandness unstranded)

Feature counts of Sample 1 (following --rna-strandness RF + -s reverse):
__no_feature        6327294
__ambiguous     2954776
__too_low_aQual     3784481
__not_aligned       688856
__alignment_not_unique      4858182

Feature counts of Sample 1 (following --rna-strandness unstranded + -s reverse):
__no_feature        6291151
__ambiguous     2911298
__too_low_aQual     4075017
__not_aligned       754400
__alignment_not_unique      16136045

Feature counts of Sample 2 (following --rna-strandness RF + -s reverse):
__no_feature        5417882
__ambiguous     1708510
__too_low_aQual     3532352
__not_aligned       564596
__alignment_not_unique      2859501

Feature counts of Sample 1 (following --rna-strandness unstranded + -s reverse):
__no_feature        5359434
__ambiguous     1676091
__too_low_aQual     3813344
__not_aligned       623122
__alignment_not_unique      2891792

These results look comparable to me across pipelines.

Thanks
Guan

If you reran the command with the correct settings, just leave it at that.  (It is not at all clear to me that strandedness rf is correct)

If you want people to tell you if you ran the commands right, you need to put down what commands you used.


Does the kinship and inbreeding coefficients depend on population frequency of an allele? I am reading Section 5.2, Kinship and Inbreeding Coefficients, of Kenneth Lange, Mathematical and Statistical Methods for Genetic Analysis. There the kinship coefficient $\Phi_{i,j}$ is defined for two relatives $i$ and $j$ as the probability that a gene selected randomly from $i$ and a gene selected randomly from the same autosomal locus of $j$ are identical by descent.
I would think $\Phi_{i,j}$ should depend on the particular population frequency of the gene or allele at that locus. However, the book does not seem to indicate that or at least does not stress that at all. Is there a dependence or not? The same question applies to the inbreeding coefficient.
More importantly, I doubt this definition is mathematically rigorous.

Here is an example to make what puzzles me clearer. It must be that the probability I am computing is not what the kinship coefficient. But this is the definition means to me. Some expert please elucidate the correct meaning of the definition.
Consider the simplest pedigree tree of one family of two parents with brothers $A$ and $B$. Allele $a$ is detected at a particular locus. Let $b$ be the other allele of the same gene. We ask for the probability of the brother $B$ having allele $a$ conditioned on $A$ having allele $a$.
In the following, to clutter symbols, we abuse the symbols by use $a$ to denote the population frequency or the unconditional probability of finding allele $a$ at the locus for the whole population as $a$. The same goes for allele $b$. $a,b\in[0,1]$ and $a+b=1$. We solve this problem by listing as follows all possible parental gene configurations, the configuration probability (considering the symmetry of paternal and maternal loci), as well as the probability finding allele $a$ for either $A$ or $B$ at the locus.
\begin{align}
\text{parent config}~~~~ & \text{Pr(config)} & \text{Pr($a$|config)} \\
aa|aa~~~~~~~~~  &  ~~a^4   & 1 \\
aa|ab~~~~~~~~~  &  4a^3b   & 1 \\
aa|bb~~~~~~~~~  &  2a^2b^2 & 1 \\
ab|ab~~~~~~~~~  &  4a^2b^2 & \frac34 \\
ab|bb~~~~~~~~~  &  4ab^3   & \frac12 \\
bb|bb~~~~~~~~~  &  ~~~~b^4 & 0 \\
\end{align}
The joint conditional probability of both $A$ and $B$ have allele $a$ at the locus is Pr($a$|config)$^2$. So the required conditional probability is
$$P(r):=\text{Pr}(B\text{ has } a|A\text{ has }a)=\frac{\sum_\text{config} \text{Pr(config)Pr(a|config)}^2}{\sum_\text{config} \text{Pr(config)Pr(a|config)}} = \frac{r^3+4r^2+2r+4r(\frac34)^2+\frac4{2^2}}{r^3+4r^2+2r+4r\frac34+\frac42}, \quad r:=\frac ab\in [0,\infty).$$
$P(r)$ depends on $r$ which depends on the population frequency of the alleles as I have claimed.
Consider $[r^3+4r^2,2r,3r, 2]$ as a positive entry vector. $P$ is a convex combination of $v:=\big[1,1,1,\frac34,\frac12\big]$. So $P$ is between the minimum and maximum of the entries of $v$, and
$$\frac12\le P\le 1.$$
We see the extrema are achieved at
$$P(r)=\begin{cases}
\frac12, \quad r=0 \\[0.7ex]
1,\quad r\to\infty
\end{cases}$$

The definition as set by Kenneth Lange's book is indeed quite vague and thus unrigorous. It is a probability not conditioned on observing any particular allele at a locus thus not on the population frequency of that allele, but conditioned on the following where we focus on one particular locus. It is rigorously defined as follows
A graph of lineage contains nodes, which we call persons, of a particular locus. Each node, or person, contains exactly two genes. We draw one and only one directed edge (parental relationship) from a gene of a parent person to the gene of a child person if the latter is inherited from the former.  This is a directed acyclic graph (DAG). Each gene can be connected by no more than $1$ edge directed to (as opposed to away from) it. The two genes in one person cannot be connected by edges. Two persons are said to be connected if there is an edge pointing from a gene of one node to a gene of the other node. Two persons cannot be connected by more than one edge. For a pair of connected persons and their directed edge between them, the edge is equally probable to connect to one gene as the other gene of a person, i.e. the probability is $\frac12$. A pair of genes (persons) are said to be connected if there exists an undirected (simple) path from one gene (person) to the other. Two thus connected genes are called identical by descent.
In this setting, the DAG with respect to the persons is given and deterministic whilst any directed edge in the DAG assigned to two persons are random with independent probability $\frac12$ with respect to the genes within each person.
For a pair of persons on a lineage DAG that is connected, the probability of one gene equally probably chosen amongst the two genes of one person being connected, or identical by descent, to one gene equally probably chosen amongst the two genes of the other person is the coefficient of kinship. It can be proven that the probability is $\displaystyle\sum_{p\in P} \frac1{2^{n(p)}}$ where $P$ is the set of all distinct (simple) paths connecting one person (not gene) to to the other and $n(p)$ is the number of persons on person connecting path $p$.
As an additional assumption, the genes that are connected have to be the same and thus assume the same allele.

With regard to the conditional probability example in the question, the coefficient of kinship is not adequate to be used to compute the desired probability. Part of the reason is the question inquires about the existence of an allele rather than picking one gene with 1/2 of the probability. We need some more fundamental probability to compute that. We will demonstrate that as follows.

We will compute the probability of the two siblings having both allele $a$ as follows.

For a pair of full siblings, they can have exactly $0, 1, 2$ identical by descent genes. Let the respective probabilities be $r_0, r_1, r_2$. Notice we are not picking genes of a person at random but talking about the existence of particular types of genes in a person. For two full siblings, some thoughts lead to $r_0=\frac14, r_1=\frac12, r_2=\frac14$.
By symmetry and without the loss of generality, we fix the connecting genes in the two siblings of the respective edges. We only need to consider the connections of the edges relative to the genes of the parents.
(1) Conditioned on there is no identical by descent genes, the sibling genes can independently assume any allele. The probability of having the configuration of $ab|ab$ (not considering gene ordering) is $2^2a^2b^2$ considering the distinct ordering of the alleles, that of $aa|ab$ is $\binom41 a^3b$, and that of $aa|aa$ is $a^4$. The total probability is thus
$$Pr_0:=\frac14(4a^2b^2+4a^3b+a^4).$$
(2) When there is exactly $1$ pair of identical by descent genes, the two identical by descent genes can only assume the same allele, while the other genes can assume alleles independently. When the identical by descent genes assume allele $a$ the probability of the two siblings both having allele $a$ is $a$. When the two identical by descent genes assume allele $b$ the desired probability is then $ba^2$. The total probability is then
$$Pr_1:=\frac12(a+ba^2).$$
(3) When there are exactly $2$ pairs of identical by descent genes, each pair can only assume the same allele but independently of each other. The only way that the two siblings not to both have allele $a$ is for the two identical by descent pairs to assume allele $b$. The total desired probability for the two siblings to both have allele $a$ is then
$$Pr_2:=\frac14(1-b^2)$$

Now we compute the probability of any one of the siblings having allele $a$. Each gene can of course assume alleles independently. So the probability is simply
$$Pr_3:=1-b^2.$$

Therefor the desired probability of both siblings having allele $a$ given that one of them has allele $a$ is
$$Pr(a):=\frac{Pr_0+Pr_1+Pr_2}{Pr_3}=\frac{\frac14a(a^3-6a^2+5a+4)}{a(2-a)}$$
since $b=1-a$. Some tedious algebra which can be carried out by, say, Mathematica shows this is identical to the expression in the question in terms of $r=\frac{a}{1-a}$.
$$Pr'(a)=\frac12\Big(2-\frac1{(2-a)^2}-a\Big).$$
It is easy to see $2Pr'(a)-2$ strictly decreases for $a\in[0,2)$ as the last two terms decreases in that interval. $Pr'(1)=0$ and thus $Pr'(a)>0,\ \forall a\in[0,1)$. $Pr(a)$ thus strictly increases for $a\in[0,1]$ and $\frac12=Pr(0)\le Pr(a)\le Pr(1)=1$.


What are the advantages and disadvantages between using KEGG or Reactome? As enrichment analysis a usual step is to infer the pathways enriched in a list of genes. However I can't find a discussion about which database is better. Two of the most popular (in my particular environment) are Reactome and KEGG (Maybe because there are tools using them in Bioconductor).
KEGG requires a subscription for ftp access, and for my research I would need to download huge amounts of KGML files I am now leaning towards Reactome

Which is the one with more genes associated to pathways ?
Which is more completely annotated ?
Is there any paper comparing them ?

One big downside of KEGG is the licensing issue. One big advantage of Reacome are various crosslinks to other databases and data.

ad 1, This depends on which pathway, they are both primary databases. Sometimes other databases that for instance combine data of primary databases have better annotation of pathways (there is an example in the review paper bellow)

ad 3, There is very extensive relatively new (2015) review on this topic focused on human pathways: Comparison of human cell signaling pathway databases—evolution, drawbacks and challenges. However I could not find there which one is more complete ...


Why is library(Seurat) causing my R session to have a fatal error? I recently tried to load library(Seurat) after installing Seurat package. It caused my R studio to have a fatal error and abort my session. I have successfully downloaded it before. I have uninstalled R and R studio and reinstalled, restarted my computer. I am new to R too, so I'm not sure what I'm doing.

I am following the directions provided on the satijalab.org website on how to install Seurat.

Here's a screenshot of the error message.

I ran install.packages("Seurat") in R. R is also crashing as soon as I run library(Seurat). Here's a screenshot right before I hit run, and when it crashes, R application just closes.

If you haven't already done so I would make a complete installation of Anaconda, particularly good with bash. The I would either install using Anaconda Navigator, or ...

# Create an R environment and activate it
conda create -n r_env r-essentials r-base python=3.7
conda activate r_env

# install Reurat
conda install -c bioconda r-seurat

# installing standard packages
conda install -c r r-randomforest --yes

If you have a Mac the instructions are here. In some systems you need to just activate r_env .

Please note if you've installed e.g. python3 with pip a conflict occurs with the Anaconda package.


Biological datasets for polythetic classification Question: What are some instances of polythetic datasets in biology? In particular, I am looking for a dataset to benchmark a machine learning algorithm optimized via episodic training.

On polythetic classification
The conventional definition of a conceptual or monothetic class is that its members must possess certain properties in common [1]. For example, in the MNIST database of handwritten digits, images can be classified into 10 distinct classes (i.e. digits) where all instances of any given class possess the same set of attributes (e.g. zeros are typically written as ellipses).
In contrast, a polythetic class encompasses members that might share a number of characteristics, but none of these attributes is essential for membership of that group or class. For example, one could group green zeros together with red ones in one polythetic class, and red zeros together with green ones in another one, as shown below.

Are there any instances of biological datasets with similar properties? I am particularly thinking of genetics data (e.g. single nucleotide polymorphisms) where perhaps a certain mutation A counteracts the effect of mutation B? Any pointers are welcome.

I may not be understanding exactly what you are talking about, but it seems to me that nearly any set of classes of interest must be polythetic. Most groups of comparable things differ in some characteristics but also have characteristics in common, and none of those characters is 100% identical with membership in the class.
Some examples off the top of my head:

Biological species / taxonomic groups (classify according to phenotypic characters)
Biological species / taxonomic groups (classify according to molecular sequence e.g. DNA)
Cell types (classify according to cell epitopes)
Cell types (classify according to transcriptional profile)
Ecosystems (classify according to species profile)
Disease state (classify according to any one of many molecular profiles)
Members of population groups (classify according to genetic variation)

etc.


Choosing Fisher's Exact or Binomial test for overrepresentation in PANTHER The PANTHER website offers a tool to obtain the GO-based overrepresentation of a gene list (the analyzed list) versus a reference gene list. After entering both these lists, the tool asks to choose between the Fisher's Exact or Binomial tests.
I guess that the convenience for one or another test would depend on the number of genes to analyze. If that is the case, which would be the ideal or advisable rank of gene numbers in both analyzed and reference lists (or their ratio) for each type of test?
(In particular I am working with a list of around 150 versus a reference list of around 2000 genes).

For a guide on using the tests in PANTHER, see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6519457/.
In answer to this specific question, the defaults are the best choice (Fisher's Exact Test, with FDR correction). Fisher's Exact test gives more accurate statistics for relatively small counts (e.g. when only a few genes in your list are associated with a particular GO term). The FDR is a better correction to use for these types of analysis in general, especially given the large overlap between classes.
The other options are there mostly to support replication of older analyses, e.g. to replicate a published analysis that used the PANTHER over-representation tool prior to 2018, when the binomial with Bonferroni correction was the default.


Extract sequences from MultiFASTA aligned file, by coordinates I am trying to extract a specific sequence from a multifasta file, from each sequence in the aligned file.
The sequences look like this, and there are 32 sequences within the multiFASTA:
>NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome
ATT AAA GGT TTA TAC CTT CCC AGG TAA CAA ACC AAC CAA CTT TCG ATC TCT TGT AGA TCT GTT CTC TAA ACG AAC TTT
AAA ATC TGT GTG GCT GTC ACT CGG CTG CAT GCT TAG TGC ACT CAC GCA GTA TAA TTA ATA ACT AAT TAC TGT CGT TGA
CAG GAC ACG AGT AAC TCG TCT ATC TTC TGC AGG CTG CTT ACG GTT TCG TCC GTG TTG CAG CCG ATC ATC AGC ACA TCT
AGG TTT CGT CCG GGT GTG ACC GAA AGG TAA GAT GGA GAG CCT TGT CCC TGG TTT CAA CGA GAA AAC ACA CGT CCA ACT

Using the following code I can pull the sequence, using coordinates, from a single FASTA file, but how would one iterate this to pull the section (by coordinates) from each sequence within the multiFASTA?
from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

record = SeqIO.read("Whole genome.fas", "fasta")

with open("Spike protein.fasta", "wt") as out:
    writer = FastaWriter(out, wrap = [])
    writer.write_header()
    writer.write_record(record[21563:25384])
    writer.write_footer()

Thanks in advance!

You can use SeqIO for the writing in addition to reading the sequences. This has the advantage that it will always write the correct fasta header, which is not always just the .id if there is whitespace, e.g. (untested)
from Bio import SeqIO

with open("Spike protein.fasta", "w") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        seq_record.seq = seq_record.seq[21563:25384]
        SeqIO.write(seq_record, out, "fasta")

BTW don't put spaces in your filenames, it can lead to issues


What does the 'suppressed' folder mean on the NCBI ftp genomes website? I wish to understand the term suppressed in folder nomenclature on the NCBI ftp genomes website?
Example: here ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/Oscillatoriales_cyanobacterium_C42_A2020_001/all_assembly_versions/
Normally, I'd see folders with names like representative or latest assembly, which are pretty straightforward. But what does suppressed refer to in this case?
Thank you for your time!

Here you have the explanation:

Suppressed genomes: Several hundred RefSeq bacterial genomes that did
not pass assembly or annotation quality validation have been
suppressed(RefSeq genome quality criteria). An FTP report file listing
the most recent suppressions will be provided soon in the RefSeq FTP
site.


How can I get or create a reference genome for Bacteria? I am a computer engineer and nowadays trying to grasp some concepts of Bioinformatics particularly, reference genomes and genomic variants. My aim is to find the effect of sequence features on variant types that occur in bacteria genomes. I chose Escherichia coli as the bacteria for my study since there seem to be more publication and data about it and its genomic sequences.
When I made a research about variants, I found out that we need to have a reference genome for the comparison with other alleles in order to find variants. And, I even read that we can use genomic sequence of one bacteria as the reference but it won't represent the cohort in an accurate manner. I got the genomic sequence data for several strains of Escherichia coli bacteria that been sequenced before. However, most of these sequences are not complete ones. I filtered out the complete DNA genomes. When I checked the reports on these genomes, I learned that these are complete genomes but not reference genomes. Only the ones which have RefSeq category: Reference Genome in their report are reference genomes for Escherichia coli.
Learning that most of the data I got for different strains of Escherichia coli were not reference genomes, I tried to seek to find reference genomes and I got this website of NCBI ref sequences for prokaryotes. This website shows only two reference sequences which match with what I have find in all sequencing data I got before. So, I want to ask two questions
First: How can I get a reference sequence for other strains of E. coli?
Second: If I cannot get a reference sequence for other strains, can I use complete genome data as the reference one? When I check the complete genome data, I cannot understand their order because there are other data for same strains that show loci of sequences but in the complete genome, I cannot see any order of the loci. Even, I cannot see data of some loci in the complete genome which blows my mind more.
I know it may not be appropriate to ask two questions in a single post, however, these two questions seem to support each other. I will appreciate any of your help.

How can I get a reference sequence for other strains of E. coli?

You can use any sequence you want as reference for a variation analysis. It just depends on the question you want to answer.
So to keep it simple, your Question is:
What are the variant region from a fresh sequenced E. coli K12 sample?
To know the differences from your known strain E. coli K12. A reason for this could be that you work with the strain K12 through out several years. During your studies you changed the conditions for the bacteria, like adding antibiotics or change the growth media. So some of the E. coli would adapt to the antibiotic and become resistant to it or they adapt better to the new media. This happens because some nucleotide changed. With each new experiment you used the best adapted colonies from the previouse one.
So you sequence your bacteria and compare it with the E. coli K12 strain. In this case this is your reference. If you use another strain, like SE15 or SE11, as reference, your result would not give you any answer for your question.
If your used data and publications gives you the information which strain were used, then use the given strain as reference.


Intersecting multiple data frames I have to intersect multiple data frame to find common elements. For visualization purpose i can use upset Plot but it doesn;t return a object containing the intersection object .I tried jvenn but at max it can take six sets to intersect .

Im able to do the interesction but i would like to know if i can do it more efficiently ,which apparently not sure how to do this is my code as of now where i read each file as data frame, covert them to list and pass them to do intersection .

NewList <- split(GO0000122, f = seq(nrow(GO0000122)))
NewList1 <- split(GO0006351,f=seq(nrow(GO0006351)))
NewList2 <- split(GO0006355,f=seq(nrow(GO0006355)))
NewList3 <- split(GO0006357,f=seq(nrow(GO0006357)))
NewList4 <- split(GO0006366,f=seq(nrow(GO0006366)))
NewList5 <- split(GO0030154,f=seq(nrow(GO0030154)))
NewList6 <- split(GO0045892,f=seq(nrow(GO0045892)))
NewList8 <- split(GO0045893,f=seq(nrow(GO0045893)))
NewList7 <- split(GO0045944,f=seq(nrow(GO0045944)))

a <- intersect_all(NewList,NewList1,NewList2,NewList3,NewList4,NewList5,NewList6,NewList7,NewList8)

b <- do.call(rbind.data.frame, a)

Can i read the files and pass them into list and do the same that would help me from doing lot of things manually

I can read the files from the directory where i have my files but not sure how to pass them to make it as list.

filenames <- gsub("\\.csv$","", list.files(pattern="\\.csv$"))
for(i in filenames){
  assign(i, read.csv(paste(i, ".csv", sep="")))
}

The object filenames is character here .It just returns me list of files.

So how to make my question bit more clear , i would like to read the list of file from the directory and then pass these to perform intersection. Any suggestion or help would be really appreciated .

Change dataframe to load all the files into a list:

wd <- "C:/all_your_dataframes_in_one_folder"

setwd(wd)

Then you get the list of files you want from your folder

files <- list.files()

Put into a list with a loop

dflist <- list()
for(i in 1:length(files)){
 dflist[[i]] <- read.csv(files[i])
}

Then use tidyverse solution

library(tidyverse)

map(dflist, ~.$GeneID) %>% #creates a list just of the column of interest
    reduce(intersect)  #applies the intersect function cumulatively to the list


Is there a Python package to convert InChi to molecular structures? I am looking for a python package that can convert InChi keys (e.g.SGNXVBOIDPPRJJ-UHFFFAOYSA-N) to molecular structure or SMILES strings. Can I do that with RDKit?

The key comes from: https://www.metabolomicsworkbench.org/databases/refmet/refmet_latest.txt

Yes, RDKit can be used, however, if you installed it with conda it will not work out of the box for inChi key fetching. You can either spend some time installing the missing bit or use something else. It is not like Mol2 support in RDKit, which shouldn't be used, it is just not worth the bother.

I would suggest PubChem's RESTFul API.

Say

 import requests
 ethanol = 'LFQSCWFLJHTTHZ-UHFFFAOYSA-N'
 r = requests.get(f'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{ethanol}/property/CanonicalSMILES/JSON').json()
 smiles = r['PropertyTable']['Properties'][0]['CanonicalSMILES']


Why does DESeq2 convert numeric columns to factor during differential expression analysis? I'm attempting to perform differential expression analysis using DESeq2 on a dataset of ~90 samples and ~20000 transcripts. I have a number of variables I'd like to test against, some of which are qualitative (factors) and some of which are quantitative (numeric). A sample of my data are below:

> head(data.phenotype.test)
# A tibble: 6 x 7
  Novogene_ID BORIS_ID     treatment queenness ovaries   elo sampletime
  <fct>       <fct>        <fct>         <dbl>   <dbl> <dbl> <fct>
1 W03BLGR     003W-blu-grn Qr12        0.0202    124.  1000  12
2 W03GR       003W-grn     Qr12        0.00210    66.8  923. 12
3 W03SL       003W-sil     Qr12        0.969     316.  1077. 12
4 W04BLOR     004W-blu-ora Qr12        0.00168    66.4  876. 12
5 W04GRSL     004W-grn-sil Qr12        0.175     210.   931. 12
6 W04YLWH     004W-yel-whi Qr12        0.984     302.  1192. 12

Now, let's say I want to test which transcripts are differentially expressed together with a numeric variable such as 'ovaries':

dds.transcript.ovaries = DESeqDataSetFromMatrix(countData = data.transcript.count.ovaries,
                                        colData = as.matrix(data.phenotype.ovaries),
                                        design = ~ ovaries)

This should work, but DESeq2 does something strange- it converts 'ovaries' to a factor, with the message "some variables in design formula are characters, converting to factors". But there are no character values in the supplied data frame, and certainly not in the column 'ovaries'! This conversion to factors is a big problem for me because, having converted a numeric column to factors, DESeq2 unsurprisingly reports that differential expression analysis cannot be performed because there are an equal number of factor levels to samples.

I have found a temporary workaround: if I reduce the data frame to just the 'ovaries' column, DESeq2 no longer converts the numeric data to factor levels and I'm able to perform differential expression analysis as normal. The problem is that I'd eventually like to perform multivariate analyses, e.g. with a design like:

~ ovaries + elo + treatment

However, doing so isn't going to be possible if DESeq2 keeps converting the numeric columns to factors, which it persistently does unless the numeric columns are supplied alone (which of course they can't be for a multivariate design).

What's going on here? Why does DESeq2 convert the numeric columns to factors, seemingly having misinterpreted those columns as characters?

NB: I have replicated this issue with different data on a different machine and the result is the same, so I don't believe this an issue with my R setup.

Edit: Here's the structure of the df:

> str(data.phenotype.test)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   64 obs. of  7 variables:
 $ Novogene_ID: Factor w/ 64 levels "W03BLGR","W03GR",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ BORIS_ID   : Factor w/ 64 levels "003W-blu-grn",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ treatment  : Factor w/ 4 levels "Qc","Qr12","Qr3",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ queenness  : num  0.02017 0.0021 0.96853 0.00168 0.17529 ...
 $ ovaries    : num  124.1 66.8 316.1 66.4 209.9 ...
 $ elo        : num  1000 923 1077 876 931 ...
 $ sampletime : Factor w/ 3 levels "0","3","12": 3 3 3 3 3 3 3 3 3 3 ...
 - attr(*, "na.action")= 'omit' Named int  1 2 3 4 5 6 7 8 9 10 ...
  ..- attr(*, "names")= chr  "1" "2" "3" "4" ...

Edited:

It's the "as.matrix" that is causing the problem.  When I use that on my sample annotation data, and make dds with that, it turns my numbers into factors


Is there a way to retrieve ENSEMBL IDs from a search query? I would like to be able to batch download FASTA files from ENSEMBL. I normally would use the API to download them from the ENSEMBL gene IDs. However, the IDs aren't ones I know, only ones I would get from searching on the ENSEMBL website and collecting the IDs manually.

Is there a way I can use the API to query search terms and collect all results that are returned?

So for example, if I use the ENSEMBL website to search "HLA" I could collect the list of IDs (ENSG00000204252 etc. etc.) manually. HLA works okay using the API to query gene IDs because there's only a few terms (A, B, C plus non-classical for class I) but I was wondering if there was a way to directly access the search query one the homepage programmatically, as I may need to do it for some messier examples.

I would download the complete geneset FASTA file from ensembl and then parse out the genes you want with a script.

Alternatively, you could do similarly with API. Get all the gene names, select the ones you want, and then pull down each sequence.


Identifying mutually **exclusive** gene sets I am interested in identifying gene pairs (or better: sets of genes) whose expression is mutually exclusive.
Ideally, both genes (or gene sets) would be widely expressed but I am also interested in cases that are very cell type specific.

Familiar examples of the former would include some of the cyclin genes, in humans e.g. cyclin E and cyfclin B, which are expressed during S-phase and M-phase, respectively. Another example would be some of the clock-genes, where Per1 and Per2 are often highly expressed in the suprachiasmatic nucleus (SCN) before midday, and whereas Bmal1 is expressed after midday.
However, such genes need not be cyclically expressed, they can simply be very specific markers for certain cell types. For example, GAD65 and GAD67 are very specific for interneurons, whereas CamKIIa is very specific for excitatory glutamatergic cells.

I have been able to find a single paper so far that systematically tries to find such genes from gene expression data, namely: Tzanis & Vlahavas (2010) Mining for Mutually Exclusive Gene Expressions. According to Google Scholar, this paper is cited by a single other paper that does not seem to be relevant. In Tzanis & Vlahavas (2010), they use a SAGE library (from 2002) of short expression sequence tags to find these strong negative associations. By today's standards, the data set is pretty small (basically, it contains 90 samples). The paper itself is hence anything but a canonical source for such negative interactions.

Questions:

Is there a curated database of mutually exclusive gene pairs or gene sets?
If not, is there a curated database that would allow me to find such genes?

Otherwise, my current plan would be to use the human or mouse cell atlas single cell data to define a signed, weighted gene expression network.
I would then search for pairs of clusters/modules, where genes have predominantly positive correlations within a cluster but predominantly negative correlations between clusters.
So essentially, signed, weighted gene correlation network analysis (WGCNA) with extra steps.

Has somebody published such an analysis already (the WGCNA literature is pretty substantial)?
Is there other (recent) related work that I am unaware off?

I am fairly new to bioinformatics, so the lack of papers and resources that I am able to find is likely mostly a reflection of my lack of ability to find the correct search term.

Addendum

Similar biostars thread from 5 months ago, also without an answer.

"Mutually exclusive" is not a precise statistical term because it could be seen as independence.

Anyway if truly mutually exclusive AND their behaviour follows periodicity, i.e. sine waves, which it might do in a cell cycle, i.e. the cells divide at periodic intervals, you can do this via asynchronicity.

If the gene expression can be described under periodicity analysis you simply sum the amplitude between two genes. When the amplitude is zero ... they are expression is antogonistic - when one is expressed the other is silent.

Amplitude simply means the height of the period.

Periodicity analysis is quite simple if you have uniform time measurements.

If you have 'missing data', i.e. you time intervals are not uniform, for example you skipped a reading then this is called missing data in periodicity analysis. Missing data periodicity analysis is complicated, but its doable.

So in this context its quite simple, you would identify the periodicity of all genes and construct a 2x2 matrix subtracting their respective amplitudes. It would be kinda juicy to do.

The above is cool if the genes are periodic, however covariance is a better solution in general terms.

What you want is negative co-variance. If two genes had a covariance of -1 ... that will mutually exclusive. As above 2x2 matrix and the co-variance across the grid. You do need to be careful that the program you are using doesn't remove negative values of the covariance.

Covariance = 1, meaning two genes are co-expressed
Covariance = 0, random distribution
Covariance = -1 mutually exclusive expression


How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome? This is a problem I have to solve frequently, and I'd be interested in knowing what other methods people use to solve the same problem.

About twice a year, I get asked to determine variants from Illumina reads, usually from either mouse or human. These are things that have good reference genomes (e.g. human genome from ensembl, mouse genome from ensembl), and are frequently tested out on various bioinformatics tools.

I have at my disposal three computers, in order of usage preference: a laptop with 4 processing threads and 8GB of memory, a desktop with 12 processing threads and 64GB of memory, and a server with 24 processing threads and 256GB of memory. Occasionally I get access to better computers when I need to do lots of processing, or to do things quickly, but these comprise my basic bioinformatics team.

I've had whole-genome projects that involved processing data from 4 individuals, up to projects with about 120 individuals, with each individual having about 10-40M paired-end reads (100-125bp each). Beyond this level, I suspect other researchers will typically have their own bioinformatics team and/or capabilities, which is why I don't get any of the larger-scale jobs that I hear about from salaried bioinformaticians.

What is the standard approach that you would use to process these reads into reference-mapped variant data (i.e. gzipped VCF) for use in downstream analysis tools (e.g. VEP)?

1. Adapter Trimming

One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it reduces the effort the mapper needs to go through to make a match (which may help some borderline reads to be mapped properly).

My preferred adapter trimmer is Trimmomatic. It has the ability to search for palindromic adapter read through, and includes an adaptive sliding window option. Trimmomatic is threaded, and seems to work reasonably fast when run through files one at a time in threaded mode. Here's an example script that I might run to trim reads from some samples. This puts trimmed FASTQ files into a trimmed subdirectory:

TRIMPATH=/data/all/programs/trimmomatic/Trimmomatic-0.36;
JARPATH="${TRIMPATH}/trimmomatic-0.36.jar";
ADAPTfile="${TRIMPATH}/adapters/TruSeq3-PE-2.fa";
CMDS=" ILLUMINACLIP:${ADAPTfile}:2:30:10:5:true LEADING:3 TRAILING:3  SLIDINGWINDOW:10:20 MINLEN:40 ";

mkdir -p trimmed;

for INFILE1 in *_1.fastq.gz; do
  base=$(basename "$INFILE1" _1.fastq.gz);
  echo ${base};
  INFILE2="${base}_2.fastq.gz";
  OUTFILE_P1="trimmed/${base}_P1.fastq.gz";
  OUTFILE_P2="trimmed/${base}_P2.fastq.gz";
  OUTFILE_U1="trimmed/${base}_U1.fastq.gz";
  OUTFILE_U2="trimmed/${base}_U2.fastq.gz";
  java -jar "${JARPATH}" PE -threads 20 -phred33 \
        "${INFILE1}"  "${INFILE2}" \
        "${OUTFILE_P1}" "${OUTFILE_U1}" "${OUTFILE_P2}"  "${OUTFILE_U2}" \
        "${CMDS}";
done;

Trimmomatic is somewhat sensitive to bad input data, and will loudly fail (i.e. stop running) if it sees quality strings of different length to sequence strings:

Exception in thread "Thread-1" java.lang.RuntimeException: Sequence and quality length don't match: 'TACATGGCCCTGAAATGACTTTCACCCAGGCAACCAGTGCCCCCTGTATAGACACATGCCTTGGGCGCTCCCCACCCTTCCTCGCGTGGCCACACCTCTGT' vs '-AAFFJFJJ7A-FJJJFJFJFJJFFA-FFAF-<A-<FFFJA-<-A7-F<<FFJJAJJJJJJJJJ--<--7-7AJ7<AAJA--J7<ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'
        at org.usadellab.trimmomatic.fastq.FastqRecord.<init>(FastqRecord.java:25)
        at org.usadellab.trimmomatic.fastq.FastqParser.parseOne(FastqParser.java:89)
        at org.usadellab.trimmomatic.fastq.FastqParser.next(FastqParser.java:179)
        at org.usadellab.trimmomatic.threading.ParserWorker.run(ParserWorker.java:42)
        at java.lang.Thread.run(Thread.java:745)

2. Read Mapping

Read mapping finds the most likely location for a read within a target genome.
There are extremely fast approximate mappers available, but these don't yet work for variant calling, and an exact mapping approach is necessary. My current preferred mapper is HISAT2. I use this instead of BWA due to the double-read issue, and because of the local variant-aware mapping. HISAT2 is made by the same computing group as Bowtie2 and Tophat2 (JHUCCB), and is their recommended tool for replacing those other programs.

  HISAT2 is a successor to both HISAT and TopHat2. We recommend that HISAT and TopHat2 users switch to HISAT2.

The HISAT2 page includes a link to a genomic index file for the human genome, including SNP variants, which I use when mapping human reads. HISAT2 is also threaded, so I run it on the files one at a time and pipe through samtools to create a sorted BAM file:

mkdir mapped
for r1 in trimmed/*_P1.fastq.gz;
  do base=$(basename "${x}" _P1.fastq.gz);
  r2="trimmed/${base}_P2.fastq.gz";
  echo ${base};
  hisat2 -p 20 -t -x /data/all/genomes/hsap/hisat_grch38_snp/genome_snp -1 \
    "${r1}" -2 "${r2}" 2>"mapped/hisat2_${y}_vs_grch38_stderr.txt" | \
    samtools sort > "mapped/hisat2_${y}_vs_grch38.bam";
done

3. Variant Calling

I currently use samtools/bcftools to do this, but would be interested in other opinions as we recently had a grant reviewer response that samtools was a dinosaur program and better approaches were available. Samtools doesn't currently work with threads for variant calling, so I save the commands to a text file and then run them through GNU Parallel. This requires the genome files to be first downloaded and indexed via samtools faidx. This produces a set of gzipped VCF files in the variants directory:

mkdir -p variants;
(for x in mapped/*.bam;
   do echo samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa "${x}" \| \
   bcftools call --ploidy GRCh38 -v -m -O z -o variants/$(basename "${x}" .bam).vcf.gz; done) > call_jobs.txt
cat call_jobs.txt | parallel -u;

As a slow, but more accurate alternative, variants can be called for all samples at the same time:

samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa mapped/*.bam | \
  bcftools call --ploidy GRCh38 -v -m -O z -o variants/hisat2_allCalled_vs_grch38.vcf.gz

If a faster processing is desired, then the variant calling can be parallelised across chromosomes using parallel and merged afterwards using bcftools norm. The implementation of this is left as an exercise to the reader.


Index a BAM file using pysam (How) can you index a BAM file using pysam?

When I tried the intuitive pysam.index I got:

import pysam
my_bam = pysam.AlignmentFile("regular_bwamem_mapping.bam", "rb")
pysam.index(my_bam)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-5-a0b3a04ecb2c> in <module>
----> 1 pysam.index(L_bam)

~/miniconda3/lib/python3.7/site-packages/pysam/utils.py in __call__(self, *args, **kwargs)
     59             args,
     60             catch_stdout=kwargs.get("catch_stdout", True),
---> 61             save_stdout=kwargs.get("save_stdout", None))
     62
     63         if kwargs.get("split_lines", False):

pysam/libcutils.pyx in pysam.libcutils._pysam_dispatch()

TypeError: object of type 'pysam.libcalignmentfile.AlignmentFile' has no len()

Oh you silly sausage, pysam.index takes a bam file name, not a python object.
import pysam
pysam.index("regular_bwamem_mapping.bam")

will index your .bam file.


A .gtf file for HG38 for HISAT2 indexing I am trying to build a HISAT2 index using Pertea et al. procedure (https://www.ncbi.nlm.nih.gov/pubmed/27560171) using their extract_splice_sites.py and extract_exons.py pieces of code. For that, I need a .gtf file of Human Genome HG38. I have obtained this file using recommendations described here http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download

The file I have downloaded is GCF_000001405.39_GRCh38.p13_genomic.gff (1,1 Gb) and when I use on it this command python3 extract_splice_sites.py GCF_000001405.39_GRCh38.p13_genomic.gff > test.ss my test.ss file is 0 Kb! I wonder if there's another .gff file that I must use and where I can get it? This is my first time working with the human genome so I am not sure if I am doing the procedure the right way.

I had the same problem and I was able to solve it with that
https://github.com/NBISweden/GAAS/blob/master/annotation/knowledge/gff_to_gtf.md


Are there any Genome in a Bottle-like resources for non-humans (especially for invertebrates)? Genome in a Bottle is an excellent resource that provides many types of DNA and RNA sequencing reads for a single individual/cell line to test genome assembly and analysis tools. For example there are Illumina WGS reads, Oxford Nanopore, PacBio, and 10X reads available for download at high coverage depth. In addition, there are pre-aligned bam files and verified variant calls.

Is there any aggregation of similar robust single-individual sequencing resources for non-human, non-mammalian metazoan species? The only way I can consider finding such resources at the moment is going through SRA projects one at a time and hoping for the best.

In addition I would like to add "Vectorbase", the curator "Dan" is a good guy (I forget his surname right now) and it's a good team. However, please read my rather somber assessment of arthropod genomics first.

We appear to be talking about lower metazoans, e.g. not protists, or yeasts (the latter I suspect are really good for your question). I don't know anything about worms.

Its important to understand that arthropods are several (soccer) divisions lower than mammals in terms of the question you are asking, except for two model organisms.

Flybase, Drosophila melanogaster is the most important species for your purposes. The second most important species within the Class Insecta will be Anopheles gambiae (malaria mosquito).

Personally I would stick to the main players because overall whilst there are a lot of insect genomes (because there are a lot of insects), their genomes remain rough because of the difficulties in assembly and requirement to physical mapping their genomes in comparison to their genome size. Also except for D. melenogaster the (incredible) tardiness in producing a draft genome. For example the idea of PacBio is in its infancy for many entomologists: and this would help them a lot in assembling their genomes correctly. I agree things have picked a lot recently, but I remain a skeptic.

Within the arthropods I would be very surprised if the crustacea, myriapods or Chelicerata (e.g. spiders) were in better condition than hexapods.

If you are interested in protists certainly Plasmodium falciparum (African malaria) has the sort of information you seek (Plasmodb).

Summary
I don't want to discourage you, insects are great, but you need to be selective if you are from a human/mammalian genetics background.


Collapse cell barcodes distribution within 1 Hamming distance I have a barcode distribution from single-cell data, e.g:

11612552 TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
...

where the first column is the number of reads associated with the cell barcode in the second column.

I would like to collapse each set of cell barcodes within 1 Hamming distance (i.e. allowing one mismatch between barcodes), keeping only the barcode with the highest read count and adding the read counts of the duplicates to it.

I was thinking to use awk and grep to filter the barcode list in an iterative procedure using regular expressions, but I think this would require too much time. I can share my code if needed.

Is there a way to do it using unix tools or already published pipelines (e.g. UMI tools)?

This is my code at the moment (it is not optimized though):

#!/bin/awk -f

{
    matching = 0;
    for (r in regexs) {
        if ($2 ~ regexs[r]) {
        matching = 1;
        counts[r] += int($1);
        break;
    }
}
if (!matching) {
    regex = "";
    for (i = 1; i <= length($2); i++) {
        l_regex = substr($2, 1, i-1) "." substr($2, i+1);
        if (length(regex) == 0)
            regex = l_regex;
        else
            regex = regex "|" l_regex;
    }
    regexs[$2] = regex;
 counts[$2] = int($1);
    }
    # so we know that 100 records have been processed
    if (!(NR % 100)) printf "." > "/dev/stderr";
}
END {
    for (c in counts) {
        print counts[c], c | "sort -k1 -rn";
    }
}

This doesn't use standard unix tools, but rather the API for UMI tools:

#!/usr/bin/env python
from umi_tools import network

cb = []  # cell barcodes
cbc = dict()  # cell barcode counts
for line in open("UMIs.txt"):
    cols = line.strip().split("\t")
    cols[1] = bytes(cols[1], "ascii")
    cb.append(cols[1])
    cbc[cols[1]] = int(cols[0])

uc = network.UMIClusterer()
CBclusters = uc(cb, cbc, 1)
cbFinal = dict()
for l in CBclusters:  # This is a list with the first element the dominant barcode
    cbFinal[l[0]] = 0
    for x in l:  # Iterate over all barcodes in a cluster
        cbFinal[l[0]] += cbc[x]
print(cbFinal)

If you give it input in UMIs.txt of the following form (note the last two lines differ from the first with a distance of 1):

11612552    TCCTGAGCACTGCATAACTCAA
9349711 GCTACGCTACTGCATAAGTCCA
8343678 CAGAGAGGCTAAGCCTGCACAT
8161950 CGTACTAGTCTCTCCGCGGCTA
8102383 TCCTGAGCGTAAGGAGCAGATC
7110298 AAGAGGCAGTAAGGAGATAATC
6626630 TAAGGCGAACTGCATAACATGT
6390489 CGAGGCTGTATCCTCTACTAGA
6210446 AGGCAGAATCTCTCCGAGTGCT
5985219 TAAGGCGATCTCTCCGTACGTG
20  TCCAGAGCACTGCATAACTCAA
2   TCCTGAGCACTGCATAACTCNA

Then it will produce:

{b'TAAGGCGATCTCTCCGTACGTG': 5985219,
 b'TAAGGCGAACTGCATAACATGT': 6626630,
 b'GCTACGCTACTGCATAAGTCCA': 9349711,
 b'CGTACTAGTCTCTCCGCGGCTA': 8161950,
 b'TCCTGAGCGTAAGGAGCAGATC': 8102383,
 b'CAGAGAGGCTAAGCCTGCACAT': 8343678,
 b'TCCTGAGCACTGCATAACTCAA': 11612574,
 b'AGGCAGAATCTCTCCGAGTGCT': 6210446,
 b'AAGAGGCAGTAAGGAGATAATC': 7110298,
 b'CGAGGCTGTATCCTCTACTAGA': 6390489}


I need some tips and suggestions for further analysis of NGS expression data (log2cpm) I am a PhD student who inherited some log2cpm data of expression data from bulk kidney tissue from a UUO(unilateral urethral obstruction) experiment that tests a new drug. The sample material consists of:

6 x Ligated Kidney (Untreated
6 x Ligated Kideny (Treated)
3 x Unligated Kidney (Untreated)
3 x Unligated Kidney (Treated)

The previous study showed that this drug was effective against fribrosis, so the aim of my study is to investigate how this drug affects inflammation and mitochondrial function in the diseased kidney.

So I have filtered my original dataset based on MitoCarta v2 DB, so investigate how the drug affected mitochondrial function. I afterwards run a PCA analysis on log2cpm data which revealed the following:

The following figure shows PC1(84.51%) and PC2(4.69%).
Red (Ligated, Untreated), purple (Ligated, Treated), green(healthy kidney, not treated), orange (Healthy kidney, treated).

I interpret this as there is a very small difference between treated and not treated ligated kidneys, however there is a small difference.

I have just started doing bioinformatics and I am very unsure which approach I should choose further, however I have though maybe do a highly variable genes analysis (HVG) see link, isolate the 500 most variable genes and run a heatmap with hieraichal clustering. I was also thinking of looking at the loadings in PC2, and take this further. The truth is that I am very unsure what is the right approach, and my research group does not have bioinformatician I can ask anymore.

I have the following questions:

What kind of analysis you guys would suggest I use to investigate
further (Hieraichal clustering, HVG etc?)
Is there any alternative preprocessing I could use before I run the PCA, to make the analysis more accurate?
Are there any packages that you can suggest?
Any other tips?

Any other tips?

I asked this same question on Biostars as well
https://www.biostars.org/p/367722/#367730

Since you indicate already having performed differential expression and using gene set enrichment the next step is to look at what changes due to drug treatment in the pathways of interest while having a decent understanding of those pathways. Things like IPA or pathview (from Bioconductor) can aid in this, but at the end of the day you need to have a decent understanding of the pathways yourself to judge what seems to be like the best targets for future work.