text:
Error generating count data using featurecounts in R
I am doing some RNA analysis and am having issues trying to generate count data. I mapped my reads to a reference genome fasta file (genbank fasta file from ncbi) using bbmap and .sam files as the output.
I am now trying to use featurecounts to generate count data but am running into issues with my gtf file. The gtf file is also from ncbi and corresponds to the alignment fasta file used for mapping. This is it's structure:
#gtf-version 2.2
#!genome-build ASM14559v1
#!genome-build-accession NCBI_Assembly:GCA_000145595.1
CP002120.1      Genbank gene    1       1362    .       +       .       gene_id "SAA6008_00001"; transcript_id ""; gbkey "Gene"; gene "dnaA"; gene_biotype "protein_coding"; locus_tag "SAA6008_00001";
CP002120.1      Genbank CDS     1       1359    .       +       0       gene_id "SAA6008_00001"; transcript_id "unassigned_transcript_1"; gbkey "CDS"; gene "dnaA"; locus_tag "SAA6008_00001"; product "chromosomal replication initiator protein DnaA"; protein_id "ADL64089.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank start_codon     1       3       .       +       0       gene_id "SAA6008_00001"; transcript_id "unassigned_transcript_1"; gbkey "CDS"; gene "dnaA"; locus_tag "SAA6008_00001"; product "chromosomal replication initiator protein DnaA"; protein_id "ADL64089.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank stop_codon      1360    1362    .       +       0       gene_id "SAA6008_00001"; transcript_id "unassigned_transcript_1"; gbkey "CDS"; gene "dnaA"; locus_tag "SAA6008_00001"; product "chromosomal replication initiator protein DnaA"; protein_id "ADL64089.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank gene    1640    2773    .       +       .       gene_id "SAA6008_00002"; transcript_id ""; gbkey "Gene"; gene "dnaN"; gene_biotype "protein_coding"; locus_tag "SAA6008_00002";
CP002120.1      Genbank CDS     1640    2770    .       +       0       gene_id "SAA6008_00002"; transcript_id "unassigned_transcript_2"; gbkey "CDS"; gene "dnaN"; locus_tag "SAA6008_00002"; product "DNA-directed DNA polymerase beta subunit"; protein_id "ADL64090.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank start_codon     1640    1642    .       +       0       gene_id "SAA6008_00002"; transcript_id "unassigned_transcript_2"; gbkey "CDS"; gene "dnaN"; locus_tag "SAA6008_00002"; product "DNA-directed DNA polymerase beta subunit"; protein_id "ADL64090.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank stop_codon      2771    2773    .       +       0       gene_id "SAA6008_00002"; transcript_id "unassigned_transcript_2"; gbkey "CDS"; gene "dnaN"; locus_tag "SAA6008_00002"; product "DNA-directed DNA polymerase beta subunit"; protein_id "ADL64090.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank gene    3163    3399    .       +       .       gene_id "SAA6008_00003"; transcript_id ""; gbkey "Gene"; gene_biotype "protein_coding"; locus_tag "SAA6008_00003";
CP002120.1      Genbank CDS     3163    3396    .       +       0       gene_id "SAA6008_00003"; transcript_id "unassigned_transcript_3"; gbkey "CDS"; locus_tag "SAA6008_00003"; product "conserved hypothetical protein"; protein_id "ADL64091.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank start_codon     3163    3165    .       +       0       gene_id "SAA6008_00003"; transcript_id "unassigned_transcript_3"; gbkey "CDS"; locus_tag "SAA6008_00003"; product "conserved hypothetical protein"; protein_id "ADL64091.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank stop_codon      3397    3399    .       +       0       gene_id "SAA6008_00003"; transcript_id "unassigned_transcript_3"; gbkey "CDS"; locus_tag "SAA6008_00003"; product "conserved hypothetical protein"; protein_id "ADL64091.1"; transl_table "11"; exon_number "1";
CP002120.1      Genbank gene    3396    4508    .       +       .       gene_id "SAA6008_00004"; transcript_id ""; gbkey "Gene"; gene "recF"; gene_biotype "protein_coding"; locus_tag "SAA6008_00004";
CP002120.1      Genbank CDS     3396    4505    .       +       0       gene_id "SAA6008_00004"; transcript_id "unassigned_transcript_4"; gbkey "CDS"; gene "recF"; locus_tag "SAA6008_00004"; note "Required for DNA replication; binds preferentially to single-stranded, linear DNA"; product "DNA repair and genetic recombination protein"; protein_id "ADL64092.1"; t$
CP002120.1      Genbank start_codon     3396    3398    .       +       0       gene_id "SAA6008_00004"; transcript_id "unassigned_transcript_4"; gbkey "CDS"; gene "recF"; locus_tag "SAA6008_00004"; note "Required for DNA replication; binds preferentially to single-stranded, linear DNA"; product "DNA repair and genetic recombination protein"; protein_id "ADL640$
CP002120.1      Genbank stop_codon      4506    4508    .       +       0       gene_id "SAA6008_00004"; transcript_id "unassigned_transcript_4"; gbkey "CDS"; gene "recF"; locus_tag "SAA6008_00004"; note "Required for DNA replication; binds preferentially to single-stranded, linear DNA"; product "DNA repair and genetic recombination protein"; protein_id "ADL640$
CP002120.1      Genbank gene    4518    6452    .       +       .       gene_id "SAA6008_00005"; transcript_id ""; gbkey "Gene"; gene "gyrB"; gene_biotype "protein_coding"; locus_tag "SAA6008_00005";
CP002120.1      Genbank CDS     4518    6449    .       +       0       gene_id "SAA6008_00005"; transcript_id "unassigned_transcript_5"; gbkey "CDS"; gene "gyrB"; locus_tag "SAA6008_00005"; note "TIGRFAM: DNA gyrase, B subunit PFAM: DNA gyrase, subunit B domain protein; ATP-binding region, ATPase domain protein domain protein; TOPRIM

When I try to run featurecounts in R using this code:
featureCounts("s1_mapped.sam", annot.ext = "GCA_000145595.1_ASM14559v1_genomic.gtf", isGTFAnnotationFile = T)

I get the following error:
The chromosome name contains unexpected characters: "CP002120.1 Staphylococcus aureus subsp. aureus str. JKD6008, complete genome" (77 chars)
featureCounts has to stop running

FATAL Error: The program has to terminate and no counting file is generated.

Error in featureCounts("s1_mapped.sam", annot.ext = "GCA_000145595.1_ASM14559v1_genomic.gtf",  :
  No counts were generated.

Do I need to format my gtf file in some way? I have also tried the code with GTF.featureType = "gene" but I get the same error.

The error specifies that the Chromosome name is not correct (I think the "," or the spaces are causing the problems). According to the GTF the chomosome is just called "CP002120.1" without the additional text. Having all this text is really quite unusual. Can you get the header of the sam file? samtools -H your_file.sam should do this or look at the first lines of the sam file in a text editor. The chromosome names start with @SQ

Also I just realized that the GTF you posted seems truncated. The line with "Genbank CDS     3396" ends with "$" as do the next two lines. This does cause the whole third line from the bottom to be quoted. If this is the case in your file and not just some copy paste error your GTF is damaged


text:
Calculating dn/ds ratio between two model organisms
I am trying to calculate the dn/ds ratio between the human genes and c.elegans. However, I am getting NA values in the output file
I used the following script
library(biomaRt)

## This gives the full list of all marts, including archives
archives <- listEnsemblArchives()

## For the c.elegans data
ensemblcelegans=useMart(host = 'jan2020.archive.ensembl.org', biomart =
'ENSEMBL_MART_ENSEMBL',
                 dataset = "celegans_gene_ensembl")

filters = listFilters(ensemblcelegans)
attributes = listAttributes(ensemblcelegans)

Celegans_GFList <- scan("c.elegans.tsv", what = character())
Celegans_GFListValues<-getBM(attributes = c('ensembl_gene_id',
'hsapiens_homolog_ensembl_gene',  'hsapiens_homolog_dn',
                                    'hsapiens_homolog_ds', 'hsapiens_homolog_orthology_type'),
                     filters = 'ensembl_gene_id',
                     values = Celegans_GFList,
                     mart = ensemblcelegans)

Celegans_GFListValues1 <- subset(Celegans_GFListValues,
Celegans_GFListValues$hsapiens_homolog_orthology_type == 'ortholog_one2one')
Celegans_GFListValues1$dnds_ratios <-
Celegans_GFListValues1$hsapiens_homolog_dn/Celegans_GFListValues1$hsapiens_homolog_ds
write.csv(Celegans_GFListValues1, file = "Celegans_GFListValues.csv")

Hi @Tom, what did you do to debug the problem? Did you check what is the result of Celegans_GFListValues1$hsapiens_homolog_dn/Celegans_GFListValues1$hsapiens_homolog_ds or what is returned by getBM?

@IIrs I am getting NA values

So, is the data available on Celegans_GFListValues1 ?

The one to one ortholog data is present but dn and ds values are coming as NA for entire input list

Ensembl discontinued our calculation of dN/dS scores for orthologue pairs in release 100 (https://apr2020.archive.ensembl.org/index.html). However, even when we did perform these calculations in Ensembl 99 and earlier, we only did this for pairs of orthologues that diverged within the last ∼100 million years, namely between pairs of mammalian, avian or tetraodontiforme genomes. More information can be found here:
https://academic.oup.com/database/article/doi/10.1093/database/bav096/2630091

Thank you so much

source: Stack Exchange

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

sort: votes


text:
Extracting Sequence Unique to a Certain Genome
I have two genomes with high similarity. But I found out that one genome has a longer sequence of about 200 kbp. I try to align the sequence with Mauve. How can I extract this unique 200 kbp?

How is this region demarcated? Between certain features? A separate contig?
How would you want the region identified for any such code to retrieve it?

I've found the answer.
I've aligned these genomes (gbk files) with Mauve Contig Mover, and then I exported the Gap File into a csv file. Then I used python pandas to analyze the region that is unique to each strand from the csv file.

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


text:
maxEE parameter from DADA2 package in R
I have just started the QC steps from the dada2 pipeline, and have failed to find a detailed explanation of what the maxEE argument entails.
I have surfed many forums, as well as the details given by the creators of the package, but they are lacking in detail. All it says is that:
After truncation, reads with higher than maxEE "expected errors" will be discarded. Expected
errors are calculated from the nominal definition of the quality score: EE = sum(10^(-Q/10))

What is Q in this instance? Is it the Quality score obtained from the plotQualityProfile function?
In the tutorial, it states that:
The standard filtering parameters are starting points, not set in stone. If you want to speed up
downstream computation, consider tightening maxEE. If too few reads are passing the filter,
consider relaxing maxEE, perhaps especially on the reverse reads (eg. MaxEE = c (2,5)), and
reducing the truncLen to remove low quality tails.

What is 2, and 5 in this instance? What does an expected error of 2, or 5, actually mean?
Thanks

The previous paragraph in the tutorial provides a link to the paper.

We’ll use standard filtering parameters: maxN=0 (DADA2 requires no Ns), truncQ=2, rm.phix=TRUE and maxEE=2. The maxEE parameter sets the maximum number of “expected errors” allowed in a read, which is a better filter than simply averaging quality scores.


text:
Correspondence of SARS-CoV-2 annotations (Nextstrain clades - Pango lineages)
Is there any annotation file providing correspondence between Nextclade and Pangolin variant nomenclature/annotations, to annotate some SARS-CoV-2 genomes from Gisaid with both these?
For the moment I could only find:

Nextclade's clade definitions
Pangolin's lineage designations

What I'm looking for is an annotation as follows:
lineage    clade
AY.43      21J
AY.4       21J
...        ...

The idea is to generate a table as follows:
genome                                 lineage    clade
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43      21J
...                                    ...        ...

The way I proceed now is the following:

Parse the Pangolin's lineage designation file to find a target genome
Look for the target genome in Gisaid (for quality check)
Parse the metadata associated to the selected genome to get the lineage for verification (only Pangolin lineage is reported)

These steps allow me to build the first part of the table:
genome                                 lineage
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43
...                                    ...

Then:

Paste the genome to Nextclade, wait for the analysis to be completed
Report the calculated clade

With this information, I can then build the table:
genome                                 lineage    clade
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43      21J
...                                    ...        ...

The idea behind this is that both lineages and clades can be used for "tagging" a genome, depending on the aim (a given clade e.g. 21J includes many lineages e.g. AY.43 and AY.4 in the example above) and that I don't enjoy manual and error-prone procedures;) Thanks for your input.

The GISAID nomenclature via Nextstrain and overall tree is presented here
(Crossposted on Nextstrain's discord forum)

The problem is that there's no one to one mapping between pango lineages and Nextstrain clades. One could generate the mapping from the Nextclade reference tree, there most tips have a pango lineage and a Nextstrain clades, and almost all pangos lineages appear in that tree. Other option is to use the open data metadata TSV from Nextstrain, same info.

A comment from a Nextclade's user satisfies my actual needs, I will post it here for future reference:

Hi - one possible source is the metadata file for UCSC’s daily-updated tree of public SARS-CoV-2 sequences:

https://hgdownload.soe.ucsc.edu/goldenPath/wuhCor1/UShER_SARS-CoV-2/public-latest.metadata.tsv.gz

It’s a big file – almost 3 million lines – because there are about 3M sequences in GenBank & COG-UK at this point and that number will probably grow. The last two columns of the file (10th & 11th) are the Nextstrain clade and Pango lineage assigned to each sequence based on its placement in the UCSC/UShER tree and the tree’s annotation of where each clade/lineage starts.

This command will make a file relating each Nextstrain clade to at least one lineage (often to a lot of lineages, especially for Delta):

gunzip -c public-latest.metadata.tsv.gz | cut -f 10,11 | sort -u > cladeToLineages.tsv

Some Pango lineages have no public sequences (only GISAID sequences, which can’t be shared at a public URL), so those lineages will be missing from that file. Some older lineages are not monophyletic on the UCSC/UShER tree, so those lineages are not annotated and will be missing. (In total, ~200 / ~1500 are missing.)

Remember not to omit the 'clarification' that in addition to Pango lineages, GISAID clades, Nextstrain clades there is also the WHO nomenclature. We could designate the Bioinformatics StackExchange monophyly, surely there is room for one more nomenclature?
Here's the WHO's response copied and uploaded from here

Just  to mention NCBI has adopted the Pango lineage nomenclature available on all of its sequences available via this link
I thought that Omicron was Pango lineage B12, or perhaps B.1.2? Oh I see it has been reclassified following further sequence data.

This information is now easily available in JSON format here: https://github.com/corneliusroemer/pango-sequences/blob/main/data/pango-consensus-sequences_summary.json. It's one big dictionary of dictionaries with Pango lineage names as outer keys and corresponding Nextstrain clades accessible via inner key:
{
 "A": {
  "lineage": "A",
  "nextstrainClade": "19B",
  "nucSubstitutions": [
   "C8782T",
   "T28144C"
  ],
  ...
 },
 "A.1": {
  "lineage": "A.1",
  "unaliased": "A.1",
  "nextstrainClade": "19B",
  ...
 },
 ...
}

It's straightforward to convert this to a TSV mapping using just minimal scripting.
Beware that not all Pango lineages pre-Omicron are contained therein. Small ones are missing currently.
Disclaimer: I created the repo.


text:
Offline Interactive Genome Map
I use the CGView Server for constructing circular genome maps. It is really a convenient tool for me to explore circular genome maps because it's interactive. I can search a certain region and it can automatically hover to that region. However, the job running on CGView Server will be deleted after 7 days. Are there any alternatives to CGView Server? Preferably something that's offline so I can save it for longer than 7 days. I've checked the website on CGView, but I think the java package for CGView is not interactive. Thanks in advance!

SnapGene Viewer: https://www.snapgene.com/features/
MacVector: https://macvector.com/
ApE: https://jorgensen.biology.utah.edu/wayned/ape/

I've been using Unipro Ugene and Snapgene Viewer (not for circular genomes, but they both have the option).
Snapgene Viewer is more user-friendly, but I think Ugene might have more tools that are free (if you want to copy-paste pieces of DNA, for instance, you can't do that in Snapgene Viewer, it would require the paid or demo version of Snapgene)


text:
What are some ways to check metagenomic bin quality?
I am new to metagenomic binning. I've used CheckM in order to estimate completeness and contamination values, and most of my bins of interest appear to have good values. My workflow was pretty straightforward (BBDuk for quality control > Megahit for assembly > vamb for binning). What else could I do / where do I go from here? (especially regarding quality control)

Some papers check 5S, 16S, 23S and ≥18 tRNA genes for completeness. For example, you can see this paper.

There are many other tools. If you want to get nice statistics and graphs I recommend Amber, there is also Anvio for manual refinement of the bins.


text:
Reading in external single cell data
I'm trying to read in an external single cell dataset from https://www.nature.com/articles/s41467-020-16164-1, but I am having trouble reading in the count matrix.
counts found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154778
This is how I've been reading in the raw counts.
cts = read.table("GSE154778_RAW.tar", sep=',', header=T, stringsAsFactors=F, check.names = FALSE)

The output is this:
> head(test)
                                                                                GSM4679532_K16733_barcodes.tsv.gz
1                                                                                xVZ\xbe\xc6pܧj\x95I\xed\030\xee\xfe\xf8ޔ\x93\xa7\xbc\xf2\023\xed{\x8e\x98)w\\\xef:ܣN\xdfrr\xc7\xf5\f\xaa\xcb\001\xfa[\002\021@|a\x91\024\xab\xc5^\017\xc0O\xbd`7\x97&m\xb9\x81\xe8R\x97{\x9f\023@\x9a\xf2\x81I\0368\x80N\027\xfa\xf8\x89fǾ6V\x98h\x9aF\023\x80\xef\x8f\xf2;f\xb1%@\xffif\004(4\xf4\037\020\xaf\x8d\032\037Y\021 1\xdd\xd9V\xc0\xde{\xf0U\017<?\xd6\xe0D\x9b\x96\x8e\xc1\xd4s\xdb\024
2                                                                                *\x9b\xee5\xbb\xb3g\xab%\xa5r\xd7?\x82\x89V\xe5\177\xd5L4g\xe7?\xa2\x96\x8b\xaa\xf0\0215D\xd4\xe1\x87~D\xd9ou\xda\xe7?~\xab
3                                                                                e\xb1\xe9ky\xb0w0Qy\xd7\xdd\xc6uaZ\xcb8\025-N\x96\xb6Y\xac\0015h\xeb\xe5\nL\xeb\xc4\177\x9cX\xbdU\016\xde\xcc\xdeAE\xa7\\w\xfcE_\001\xbb\026\xba\xe7\xf7ܳS'\x96\u038bI\xfe\x9a\xef\xfdm\xa2\xbf_瘺a\030\037\xec\xa3\032\x97\xa2\xea\x9e\xf3\xf40\x8d\xfb#\n\xa9\xb8vE\x94\x9b:\xef\xcb\xe5\xd5\xc69n\xe1\xbd\xcab\025\xc7\xfb|\fT\xb5\xfa\xea4lK\xb2\037)\xdfL\xfc\x8f\xf5yz\xa6D}\001/\xaak?̑\xe9\003\xf1\031\x94\xed~\xfdӳ\xc7B]\xca/\\\xfb\xa3%n
4 \xef\xaa\xc4B\xe0]\xe5\xac6\x9e\xea\xbbe\xf9\x8a\xecch/nV5\xd1p_O\xb6N\034\xa9\xf3x\034\x98&\025\xd75\xa7b\021\xb3v\xd6P\xf7D4~u\xb3\xab(\xf9\x9d\xb29^,\xd8\xf3v`B)\xc3c|>'\xc36Z\016\xb4\0257x8\u07b5fb\xd1(\xb9\xbe\xe7WT.x\xaa5\xe4\xc1\xd6\xe6ny\f֟\fYST.\\\xd0\xf2h\xf6\x8d\xc5Q\x90o\xbc\xadG\004\xb1\x9el\xcb\xd8\037\xef\xd4\024Y\xf1Q.\xe3\x81FGZ\016\026t\xa4\xcaY\xf3\xb7\x95\017,\au\xe8\037\xcf\a\xb3\ab\001\x8d.\aճ\x8f\xe6\xc4h\xe5n8\xf7Sdg\xc5Z\xef\xaa\xcb!\xd6\xd0ӳ\x92\xa2\xc4Wݥ?\xcc\xc4\016\031\023\xaa\xe7\xb6-H^X\005\xb1ӡ\x9b\xe6\xfb\xf0\x9b\x9d\027q\xcbq\x9aEm%\036\xb3\xe9$f\xc8\xd6\005\xdb\xee\x8b\xec\xf5\xb9\xef'\xaa~\xad\xb3\xf9\x8fg\xb3\025\xe5h\xc9\xe9_\xee7v\x87\xa3\xf5\177,ǌ\x9c\016Z\f\xe6\xf353\027\0367=\xc0\002{\x89\037,\017\xd3Œ\030\034D/\xeeb\x8ct\033\xe5UD\xedϧ\xb9\xf1\xbe\xe2\034L\xa7o\xea\xff\xb9\xecŭ\x9b\036Õ\xdbM+>Su$\\\xe0\xda2\xdbl?\xa2\xcc\v\xf4rS\xe4\027\x8aEdR9f\xb9N\022\xb4И@,e\xfe<\x95gN\xb4Ć\x9b\xc7\xeb\032\xc7P~\xa2\034\xd8]\x98\032\te\xff\xd8!vr\xb0\xa6\x92K\xf0\xfd`*\xb3\xe8\aq\x90\x9e'\xe6FI\x84\xa3<B\xeddE\xefpg\001V[\xbc\xb7u\xb04\003\xe3K\xb6h\xc0ild`\x92\025W\xe5#J\177\xa6\x95\xabZ.\x96\xec\016w\x93d÷+]Di\xbaJ\xa2\003\xf5\021m\xa4\x80J\022\f\027QqE\021r\x94\xf0\xbd\037\xd6ye\x970Fm\xa1ˎO\035\xfe\xe5\xb8\xf2/\x89\x9bX$\xa9b\026\037B~\xe2\023k\033\xbaB\xbcɪ\xa7\xe4f\xf2\027\x8f\xa2t\xa1\xee\xbd\xec\x80\xe5;\037\xbc\xfeC\x94\xed\x91UQ\x81~\xbe^N<\003&\xceDJE\xab\031\xa7\xf9\x9a\xa2\xf2%\024}\x9e.g\x9a0J&lj\xd2\032\xc55~\022\xa3ޭ\xd3|\xd1,k#\017ֶ\xfbg\026CL\016\x9bx\023Y\x89\xffI\x85\x8c\x89>U=\024|í\xb4\xf2\006\xfa#ʎ\032v\xc6Hps\034\fI\xef`E\xb1\xec\xcb\xf3\xc6\xf2E\xe5ӹ\xee1σc\x8b\x89\xe5\xf2\xd0\xe3\xf2\xeb\xac\xe5\xfb,&N\002\x9b\xe14i\xc5\024y\xc8^Ď\xe73\n\006\022\x94TM\xc5p\xbf(\xb8w~ES\x8d&N\xafr\xefL\f\xbf\x816~$\x99\021A\xfe7m\tƘhX\xe0\xa2\xc2>&\xe5s\xa0\x99\t\004\x88\u07bd\xa34\xf6q\xe4\xf1ٕ\x9b\xea\003\xd7frq(~\xe7\024\035\xcf
5                                                                                \x91\xb7U
6                                                                                \x9f-\xcf\025\xa7\xeb6<\x89\xf55\xc6+\xb9K\xb5\xbaa-[.\xda \xca]<M\xea\xbd%

Any suggestions?
Update:
I unzipped the file in HPC but now I'm having trouble reading in the counts as it's all the counts per sample.
Data_Processing.nb.html            GSM4679535_T3_genes.tsv.gz     GSM4679540_T9_barcodes.tsv.gz      GSM4679544_Y00014_matrix.mtx.gz
Data_Processing.Rmd                GSM4679535_T3_matrix.mtx.gz    GSM4679540_T9_genes.tsv.gz         GSM4679545_Y00016_barcodes.tsv.gz
GSE154778_dgeMtx.csv.gz            GSM4679536_T4_barcodes.tsv.gz  GSM4679540_T9_matrix.mtx.gz        GSM4679545_Y00016_genes.tsv.gz
GSE154778_RAW.tar                  GSM4679536_T4_genes.tsv.gz     GSM4679541_T10_barcodes.tsv.gz     GSM4679545_Y00016_matrix.mtx.gz
GSM4679532_K16733_barcodes.tsv.gz  GSM4679536_T4_matrix.mtx.gz    GSM4679541_T10_genes.tsv.gz        GSM4679546_Y00019_barcodes.tsv.gz
GSM4679532_K16733_features.tsv.gz  GSM4679537_T5_barcodes.tsv.gz  GSM4679541_T10_matrix.mtx.gz       GSM4679546_Y00019_genes.tsv.gz
GSM4679532_K16733_matrix.mtx.gz    GSM4679537_T5_genes.tsv.gz     GSM4679542_Y00008_barcodes.tsv.gz  GSM4679546_Y00019_matrix.mtx.gz
GSM4679533_Y00006_barcodes.tsv.gz  GSM4679537_T5_matrix.mtx.gz    GSM4679542_Y00008_genes.tsv.gz     GSM4679547_Y00027_barcodes.tsv.gz
GSM4679533_Y00006_genes.tsv.gz     GSM4679538_T6_barcodes.tsv.gz  GSM4679542_Y00008_matrix.mtx.gz    GSM4679547_Y00027_genes.tsv.gz
GSM4679533_Y00006_matrix.mtx.gz    GSM4679538_T6_genes.tsv.gz     GSM4679543_Y00013_barcodes.tsv.gz  GSM4679547_Y00027_matrix.mtx.gz
GSM4679534_T2_barcodes.tsv.gz      GSM4679538_T6_matrix.mtx.gz    GSM4679543_Y00013_genes.tsv.gz     Lin_et_al_GMedicine_metadata.csv
GSM4679534_T2_genes.tsv.gz         GSM4679539_T8_barcodes.tsv.gz  GSM4679543_Y00013_matrix.mtx.gz
GSM4679534_T2_matrix.mtx.gz        GSM4679539_T8_genes.tsv.gz     GSM4679544_Y00014_barcodes.tsv.gz
GSM4679535_T3_barcodes.tsv.gz      GSM4679539_T8_matrix.mtx.gz    GSM4679544_Y00014_genes.tsv.gz

I've created this loop to read through all the samples.
counts_list <- list()
for (sample in samples){
    print(sample)
    filename_header <- sample
    barcode_file <- file.path(wd3, paste(filename_header, "_barcodes.tsv", sep=""))
    gene_file <- file.path(wd3, paste(filename_header, "_genes.tsv", sep=""))
    matrix_file <- file.path(wd3, paste(filename_header, "_matrix.mtx", sep=""))
    counts <- Matrix::readMM(matrix_file)
    barcodes <- read.csv(barcode_file, sep="\t", header=F)
    barcodes <- barcodes$V1
    colnames(counts) <- barcodes
    genes <- read.csv(gene_file, sep="\t", header=F)
    genes <- genes[,-3]
    colnames(genes) <- c("ENSG", "Name")
    rownames(counts) <- genes$ENSG
    genes <- genes[!duplicated(genes$Name),]
    counts <- counts[genes$ENSG,]
    rownames(counts) <- genes$Name
    counts_list[[sample]] <- counts
}

Receive this error when I run:
counts <- do.call(cbind, counts_list)
Error: Matrices must have same number of rows in cbind2(argl[[i]], r)

I added:
genes <- genes[,-3]

to ensure that all the gene file has only 2 columns as I noticed that one sample had 3 variables but I'm still erroring.

Tar is an archive format so as you figured out by now you have to unpack it.
The files you now see are as follows:

mtx.gz are compressed count matrices which you load into R using the Matrix::readMM() function
genes.tsv.gz are the associated rowwise data so the "genes" which you can load with read.delim and then use rownames() to assign them to the matrix generated by readMM()
barcodes.tsv.gz are the associated column data so the samples or "cells", use colnames() to assign them to the matrix

I would use list.files to create a vector of files and then either a loop or lapply to iterate over them for the loading. Does that make sense to you?
Edit after OP received this dimensions error: This then means that these files come from different experiments and have different upstream processing. THere is nothing you can do about it but figuring out which files come from which toplevel experiment and then only combine those of one experiment, respectively. The procedure towards reading data is as I lined out above.

I added a loop in order to read the files but receive the above error (updated my question).

@mmpp This then means that these files come from different experiments and have different upstream processing. THere is nothing you can do about it but figuring out which files come from which toplevel experiment and then only combine those of one experiment, respectively. The procedure towards reading data is as I lined out above.

It looks like it is 10X data and I strongly suggest a dedicated single cell analysis package such as Seurat. You will most likely end up using one of these anyway for data analysis.
Read10X() can be a good start. I don't remember whether it requires dedicated folders per sample though. Even if this is the case, you can create individual sample folders with a simple bash script, can be done within R as well.


text:
Plink error: --check-sex/--impute-sex requires at least one polymorphic X chromosome locus
I am running an analysis in plink.
This is my first attempt at a QC and I keep encountering the same message when I try to go through with my sex check. I've done the SNP missingness step, it's this step that keeps stumping me. I'm assuming it has something to do with the chromosome naming but I'm a bit lost as to how I should tackle that problem. For reference, I attempted running the --check-sex function on the files (which are in bim, bed, fam) however, I'm met with this error.
431911 variants loaded from .bim file.
13635 people (7170 males, 6465 females) loaded from .fam.
13635 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 13635 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.997824.
431911 variants and 13635 people pass filters and QC.
Among remaining phenotypes, 4457 are cases and 9178 are controls.
**Error: --check-sex/--impute-sex requires at least one polymorphic X chromosome
locus.**

Does anybody know what I should do to fix this?
This is the command I used: plink --bfile ppmissingness2 --check-sex

Without X chromosome variation, sex cannot be reliably determined from input genotypes.
The error is telling you that there are no polymorphic X chromosome loci in the provided data. Either your input files don't include any X chromosome information (this could happen due to the increased statistical challenges involved in X chromosome data), or all the X chromosome "variants" are invariant / homozygous.


text:
How to chose values for QC for a GWAS?
I am running a GWAS on a large dataset. I am using the H3aBioNet QC workflow (https://github.com/h3abionet/h3agwas/blob/master/qc/README.md).
I have run the QC on initial parameters:

MAF < 0.01
GENO > 0.05
HWE 0.0000001

My question is now that I have the output, do I need to change these parameters to a better estimate?
I have the following graphs:
Figure 1 SNP missingness: For each level of missingness specied on the x axis, the corresponding y-value shows the proportion of SNPs which have missingness less than this.
Figure 2 Missingness per indvididual: For each level of missingness specied on the x axis, the corresponding y-value shows the proportion of individuals which have missingness less than this.
Figure 3 Cumulative frequency of SNPs. For a frequency shown on the x-axis, the corresponding y-value shows the proportion of SNPs with frequency at least this frequency; that is, it shows the proportion of SNPs which will remain if the MAF lter of this x-value is chosen.
Figure 4 HWE distribution. For an HWE-value shown on the x-axis, the corresponding y-value shows the proportion of SNPs with HWE p-value at least this frequency; that is, it shows the proportion of SNPs which will be removed if the HWE lter of this x-value is chosen.
Figure 5 QQ plot for Hardy-Weinberg scores.
Please let me know which one is the best to attach to the question.

How many SNPs do you have before and after applying those filters? If you’ve ended up removing most of them, then I’d suggest maybe rethinking what values you choose, otherwise I don’t think it necessary. Is there a reason why you thought you should?

@user438383 I think its due to the graphs that were produced. I dont know how to interpret them and I cannot add them all to the question here. There are 295218 SNPs, 80370 SNPs were left over, too stringent?


text:
Nanopore variant calling
So far I haven't done any variant calling as such. Nanopore I have used for 16s microbiome data.
Now My question/doubt so how do I proceed for nano-pore virus sequencing data
Steps:

I get fast5 files
The i perform base-calling using guppy
Get fastq files

Now after doing the above is it possible to use variant calling that is used for Ilumina once I get the fastq files?

Tools for illumina are not really suitable for long read sequencing. For nanopore data I would suggest LongShot, Clair3 or DeepVariant-PEPPER.

Sure it is possible, but I guess it won't be more reliable than a method specifically designed for Nanopore. But if you use Nanopore I guess accuracy is not your main concern.


text:
Understand this Kmer plot from Merqury?
The attached figure is generated based on Illumina reads from multiple individuals compared to genome assembly. Looks like there are a lot of kmers are reads only (grey colored). Also, a blue peak (2x kmers), does it represents high heterozygosity? GenomeScope gives 1.5% of heterozygosity result. But it seems the blue peak normally should be near the 70x on the x-axis. Thanks.

Assuming a diploid genome, usually, we aim to build a monoploid reference, representing a single haplotype.
Practically speaking, we need to include only one of the two k-mers that are representing heterozygous sites, however, we would desire all the homozygous k-mers to be included exactly once.
Looking at your plot, I would guess 1n coverage is somewhere around 18x, and 2n about 36x. If that's right, your plot is suggesting that your assembly contains lots of uncollapsed haplotypes, which is quite understandable given very high heterozygosity you estimated with genomescope. This is also supported by many homozygous k-mers are there twice (the blue bump). These are probably regions in between heterozygous sites, or flanking them.
Anyway, I would recommend something like purge haplotigs to collapse haplotypes into single monoploid reference and redo the plot. The blue bump should disappear and the black and red distributions should have about the same height around 18x, while nearly all the 36x k-mers should be red (once).
Also, if the interpretation of the plot still would not be clear, you can check this part of a lecture by Lucía Campos-Dominguez on the topic (it's a different software called KAT she talks about, but it is using the same principle as Mercury).

Hi Kamil, thank you for your informative reply. You are exactly right about the 1n = 18, although not sure how you get that so accurate ;-). I do have another Merqury kmer plot with a very low blue peak near the 4n location rather than the 2n location, while the red peak at the same 2n. And this happens for another species genome assembly. Is that normal? Because that makes the blue peak at 2n look weired.

I will post another question. Thank you!


text:
Pymol python script for selecting a residue's sidechain and calculating its center of mass
I trying to write a pythonic PyMol script that loads a PDB, cmd.iterates through the PDB and stores the residues (three letter code) into an array.  I then want to loop through this array and calculate the center of masses (COM) of each aa's sidechains (except gly) and add it into a second array.  I am able to calculate the COM of a residue as a whole, i.e. including backbone atoms.  But what I want to calculate is the COM of the sidechain.  Please point me in the right direction.
path = 'Path to PDB dir'
for filename in glob.glob(os.path.join(path, '*.PDB')):
    cmd.load(filename)
    stored_residues = []
    cmd.iterate('name ca', 'stored_residues.append(resn)')

    stored_residues_COM = []

    for i in stored_residues:
         stored_residues_sidechain_COM.append(centerOfMass('resn %s' % i))

I am guessing you used the command pymol.cmd.centerofmass? If so the selector is the ticket to doing it by selecting everything of the residue bar the backbone atoms, e.g. f'resi {r} not name N+CA+C+O+H'
Parenthetically, the cmd.iterate command is really ugly and counter-intuitive —it's heavily used because it is older. The command get_model(selection) in pymol2 module is so much nicer.
for atom in pymol.cmd.get_model('name CA').atom:
     atom.resi #note it is type str due to insertion codes!
     atom.resn
     atom.coords # List[float]

Second parenthesis. A and not element H is not needed or beneficial to pymol.cmd.centerofmass as this takes account of mass and occupancy. There is no need to select only heavy atoms as a proton weighs nothing. The barycentre and centroid are the same in geometry where there is no mass to the points, but here mass is present.

That is good advice.  Thank you.  I am also looking at Thomas Holder's code (https://pymolwiki.org/index.php/Sidechaincenters) as a comparative approach.

@gwot Upvoting and checking 'answer resolved' is a good way to say thanks

Hi Matteo, sorry for the breach of etiquette.  I have accepted you helpful answer, but do not yet have the rep to upvote.


text:
How to extract a subset of individuals from a .bgen file?
I have 10,000 controls and 3000 cases and the .bgen file has genotype of 200,000 individuals. I would like to extract only those individuals from the .bgen file that are either cases or controls using plink for downstream analyses. How can i do that?

You can use qctool for this. If you create a list of IDs in list.txt (one ID per line), you can filter the .bgen file(s) by running
 qctool -g yourdata.bgen -s yourdata.sample -incl-samples list.txt -og yourdata_filtered.bgen

what does .sample file contain?

It's a type of file in the Oxford format that usually contains IDs, missingness, sex, and other covariates (similar to a plink's .fam file). More info here: https://www.well.ox.ac.uk/~gav/qctool/documentation/sample_file_formats.html

so in my .sample file do i only need the ids that need to be extracted or all the ids present in the .bgen file?

the order of the ids must be the same as in the genotype file, how can i know the order in the genotype files?

You should have a .sample file together with the BGEN data already. If you don't, it might be worth checking with who provided the data in the first place, since you can technically obtain one from the BGEN file but it's a laborious process.

It might not have the extension .sample, but it should be a plain text file with multiple columns structured as above


text:
About workflow input section in DRAGEN-GATK v1 whole genome analysis script
I am trying to understand the WDL script for the DRAGEN-GATK whole genome analysis and I am faced with this input section inside the workflow definition which I really don't get:
# WORKFLOW DEFINITION
workflow WholeGenomeGermlineSingleSample {

  String pipeline_version = "3.0.0"

      input {
        SampleAndUnmappedBams sample_and_unmapped_bams
        DNASeqSingleSampleReferences references
        DragmapReference? dragmap_reference
        VariantCallingScatterSettings scatter_settings
        PapiSettings papi_settings

        File? fingerprint_genotypes_file
        File? fingerprint_genotypes_index

        File wgs_coverage_interval_list

        Boolean provide_bam_output = false
        Boolean use_gatk3_haplotype_caller = true

        Boolean dragen_functional_equivalence_mode = false
        Boolean dragen_maximum_quality_mode = false

        Boolean run_dragen_mode_variant_calling = false
        Boolean use_spanning_event_genotyping = true
        Boolean unmap_contaminant_reads = true
        Boolean perform_bqsr = true
        Boolean use_bwa_mem = true
        Boolean allow_empty_ref_alt = false
        Boolean use_dragen_hard_filtering = false
      }

What do the first 5 lines under input mean? Are they Strings? Files? Or what exactly?
Thanks in advance for your help
Nawar

GATK is Broad's VCF
DRAGEN is Dynamic Read Analysis for GENomics to facilitates hardware-acceleration of analytics and processes including conversion, compression, mapping, alignment, sorting, duplicate marking, haplotype variant calling and joint genotyping.
https://gatk.broadinstitute.org/hc/en-us/articles/360045944831-DRAGEN-GATK


text:
Searching for adapter sequences in FASTQ files - metgenomics
I have recently received some metagenomic data from 16S rrna sequencing. The sequencing company claim to have removed primers, however not adapter sequences. Please note that the files have been demultiplexed.
In the platform, we have forward and reverse strand fastq files. For each fastq file, there is information regarding the adapter used, for each strand.
However, when I run the following command:
cat myFastQ_forward.fastq | grep 'adapter_forward', it returns nothing. My question is, shouldn't the adapters be identified successfully simply by looking for an exact match in the fastq files? If not, why?
Thanks,

Just run it through fastqc and you will know whether adapters are present or not. Exact grep makes little sense as sequencing errors and truncation (when read length is shorter than the full-length adapter) do not allow naive greping.

It depends on the read layout they used, but typically when I hear anyone talk about removing adapters in paired-end sequencing, it's the occurrences of adapter sequences at the ends of each read from the complementary strand's adapter.  For example see the figures on this IDT page about 16S sequencing.  If a given sequence is short enough/number of sequencing cycles high enough, you'll end up seeing "past the ends" on each read.
For example, with a reverse read reverse-complemented and aligned with a forward read, you could maybe see something like this off the sequencer:
                       |-->
Forward:               BarcodefwdAmpliconBarcoderevAdapterrev...
Reverse:  ...AdapterfwdBarcodefwdAmpliconBarcoderev
                                               <--|

And then if the company already removed barcodes and primers at the start of the reads, maybe they gave you:
                                 |-->
Forward:                         AmpliconBarcoderevAdapterrev...
Reverse:  ...AdapterfwdBarcodefwdAmplicon
                                     <--|

So the material you'd want to trim is on the ends, and reverse-complemented.  How much you see depends on the lengths involved so maybe you'll see partial adapters or none at all, and there may be sequencing error, too.  See cutadapt's documentation about trimming 3' adapters and paired-end reads for more on that.


text:
How to retrieve fasta sequence after local blast?
I have created a Blast database using a reference genome. Then, I have performed a local blast search in command line using a gene of interest. I have obtained some hits with the usual Blasting information. Now, I want to extract the exact matching sequence from the blast search with their corresponding start and stop position. I have a .gff file for the genes as well for the reference genome. Can you give me any suggestions on how to do that? Thanks a lot.

Please [edit] your question and add an example of the relevant files so we can understand what you need. Do you have multiple HSPs you want to combine? Do you need to model splicing? If you already have the gff of the target sequence, why were you blasting?

FWIW, In the Sequenceserver BLAST software, you can just click on "download alignment" to extract the aligning region.

My usual approach for extracting sequence from a reference genome is to use samtools faidx:
# create genome index; slow; only do once per reference genome
samtools faidx reference_genome.fa
# extract sequence for a specific match
samtools faidx reference_genome.fa chromosome_name1:posStart1-posEnd1 > match1.fa
samtools faidx reference_genome.fa chromosome_name2:posStart2-posEnd2 > match2.fa

You haven't mentioned how many BLAST results you have, or what format they're in, so I'm not sure if it's practical to do this manually, or if it needs to be automated in some way.


text:
Chou-Fasman Method Overlapping Regions Problem - Secondary Structure Prediction
According to the Chou-Fasman method, when there is an overlapping region on the amino acid sequence, we compare the scores of alpha-helix and beta-sheet and chose the highest one.
Let's say I have this amino acid sequence: ABC.
As a result of computations, AB is identified as alpha-helix, and BC is identified as beta-sheet.
This means, I have an overlapping region, which is B in this case.
The algorithm says that we need to compare the scores of the overlapping region, and select the secondary structure with the higher score, let's say alpha-helix in this case.
Would the resulting secondary structure be $\alpha$$\alpha$$\beta$ or $\alpha$$\alpha$_?
In other words, do we cancel out the secondary structure identifications when we encounter overlapping regions?
Thanks in advance.

Sure, an isolated residue has a ψ and φ angles, so you could technically have an isolated residue in a given SS conformation, but that does not mean much if its neighbours' backbones aren't in that same conformation. Also, sure, in Rosetta design algorithms residues are thrown into ABEGO bins based on ψ and φ angles alone and as a result one can get isolated α-helix bin residues (in what is a loop): but that's a design algorithm that will want to correct this.
In your particular example, IIRC that algorithm uses a window of 4 or 6 residues for SS prediction. So it will surely be (αα)αα_.
But I am pretty certain you cannot have a ααααββββ switch without a loop or turn in nature —I might even hazard a guess that for design it is not Koga-compliant even if not specifically addressed. And with the Chou-Fasman algorithm you'd need to have the windows spanning the switch score above the threshold in the correct way. I am sure there's a combination that allows it as after all we're talking of an ancient 70s prediction method based on probabilities alone...


text:
Methylation data: beta values are normally distributed
I am trying to incorporate Methylation data (Illumina 450K) into my pipeline. I was provided with two versions of the same data: one (dataset $a$) normalized with ComBat and the other (dataset $b$) normalized using the residuals method (i.e., taking the residuals of the models fitted on the covariates). Looking at the distribution of the beta values, I noticed that:

The values for dataset $a$ have a typical bimodal distribution with peaks at around $0$ and $1$
The values for dataset $b$ are roughly normally distributed with the mean at around $0.5$.

I am not that familiar with Methylation data, so maybe this is expected. Is it? I did not find any paper showing similar results.

$\beta$ values are the ratio of methylated to unmethylated probe intensities.
This means that $\beta \in [0, 1]$, so it can't be normally distributed, by definition.
Its distribution may have sections that look like bell curves, but it isn't normally distributed.
You can take a look at this paper by Du et al. [1] for details about the statistics of $\beta$ values.
That said, the distribution of all $\beta$ values from an array or bisulfite sequencing (or some other DNA methylation assay) should be bimodal with modes at 0 and 1.
Nearly all CpGs in the human genome are methylated or unmethylated.
There may be some hemi-methylated sites, but for the 450K array you mentioned, there should be few of these.
It sounds like you might be looking at some residual value, or some other value from the analysis in the paper for dataset $b$ other than $\beta$.
Could you be looking at $M$ values, or something else?
[1]: Du et al., BMC Bioinformatics, 2010

It's not clear from your question what the statistical design of the experiment is, so I can't judge whether or not a distribution around 0.5 would be expected from the residuals after fitting.
Given that dataset a looks like what you expect, I'd be inclined to use that (and perhaps ask the person who gave you this data why dataset b looks so different, and why it ended up with a 0.5 mean).


text:
Identifying Novel Species with ANI Taxonomy
I have sequenced a bacterial genome and after I checked the ANI scores, the scores against similar bacterias were 79 and lower. How can I confirm if this bacteria is a novel species/strain? Is there any threshold value for ANI score? I'm pretty sure if the ANI score is above 95, it indicates the same species.
Should I double-check this result by using other methods such as GGDC?
Thanks in advance!

Short answer:
Yes.
Long answer:
I suggest 4 methods for checking the genome.

Use minimap2 or a similar aligner to perform alignments with the most similar genomes, then generate a dot plot or other visualization of the regions that align to figure out what is going on. Here is one tool, I haven't used it, here is another. I generally use minimap2 to very quickly align, and then a quick plot in R using the PAF file.
Use k-mer hashing methods such as mash or sourmash to quickly compare to a large database. Large precomputed databases are available for download on the linked mash site. You can also use it for one-by-one comparisons.
Use QUAST to estimate more informative summary statistics for comparisons to the more similar genomes. QUAST has a rich set of possible visualizations that you can use. It also tries to identify the most closely related public genome using rDNA sequences, though that has been hit or miss for me.
Use a quality assessment tool such as CheckM or BUSCO to estimate the quality of your genome. THese also have utilities for taxonomic placement of your genome in case something bad happened (e.g. sample swap).


text:
Extracting and analyzing structure data of multiple proteins from PDB using Python
I am new to programming and I am given the task to extract structural data of soluble proteins from PDB and then shortlist the soluble enzymes from the extracted data set for further classification based shortlisting of soluble non membrane enzymes. Since I do not have much expertise in programming I am finding it difficult to figure out how to solve the problem. Any help would be really appreciated.

It sounds very much that you are interested in enzyme properties that aren't the cartesian coordinates of the atoms. Sure, you can calculate the Gibbs free energy of folding of a structure, a proxy for stability, but converting this to melting temperature without a reference is crude and at that point you may as well use the ProTherm database. One can also do analyses of the Poisson-Boltzmann surface electrostatic charges etc. but for solubility simply calculating the isoelectric point from a primary sequence will be easier and will circumvent the need for a structure...

@MatteoFerla Thank you for your response. Is there a way by which I can import the pdb files of proteins using tags like Soluble, Membrane/Non-Membrane ? The initial Import of soluble proteins from the database is an issue that i am facing. Once i am able to import all the proteins, then I plan to use Binary SVM to further classify then into enzymes and non enzymes

The classification of enzyme vs. non-enzyme should definitely not be done with a classifier on extracted structural properties. Parsing the UniProt entry for the protein has all you need. Take an enzyme: https://www.uniprot.org/uniprot/P06721 there's the field in the XML (sorry no JSON) that gives you the EC number and another the subcellular locationation. It also gives what PFam domains are there. Alternatively one extract the GO terms. There's no pI, but there's a ProtParam mimicking function in BioPython for that.


text:
10X low rate of correct barcodes was observed for the candidate chemistry choices for the input
I am testing a 10x fastq dataset , but cellranger count complains "An extremely low rate of correct barcodes was observed for all the candidate chemistry choices for the input", I have tried all the chemistry recommendation codes on the 10x genomics website, and I still get the same result, how can I solve this problem, thanks.

It is difficult to provide an answer when the question does not include what the input is, what command was run, ... I can only suggest to use the --chemistry=auto option of cellranger count.

What does the quality look like on those reads?  If this is 10XGenomics library, they should be able to help your troubleshoot.

Could you please give a few more details about the nature of this project? Are you sure this is a single-cell dataset (I'm aware that 10X do other non-single-cell things as well)? Showing the first few lines of header from the fastq files can be very helpful to work out the structure of the experiment.

how were your fastqs generated? and how are they named?

Thanks @ciwei for providing your answer. Could you kindly check the "accept" box for your answer? This does help site statistics

I have solve this problem, cause the combined barcode is differernt from the common 10x barcode. when mapping the genome, need to add-- chemistry ARC-v1 when cellranger count.

Excuse me, how do you determine the required chemistry?

ARC-v1 means this was done using a multiomic kit. From cellranger count --help "To analyze the GEX portion of multiome data, chemistry must be set to 'ARC-v1'; 'ARC-v1' chemistry cannot be autodetected"

Answer from @maximilian-press, converted from comments:

I can't comment on 10X data specifically, however, I suggest writing a
simple script or command that counts the frequency of each barcode,
and reports back the N (10? 20? 100?) most frequent ones, which will
allow you to investigate in more detail. This is what I have done in
the case of suspected bardode swaps or reverse complement issues.


text:
HGVS python package : How to parse complex insertions?
In HGVS package, how to create more complex variants ?
I just wonder how can I parse more complex insertions such as those :
NM_004006.2:c.940_941ins[903_940inv;851_885inv]
NC_000006.11:g.8897754_8897755ins[(543);8897743_8897754]
NC_000006.11:g.10791926_10791927ins[NC_000004.11:g.106370094_106370420;A[26]]
(All taken from the varnomen website)
When I tried I got parsing errors (mainly because of the brackets that represent the inserted sequences)
Whole traceback :
In [6]: var_g = hp.parse_hgvs_variant('NM_004006.2:c.940_941ins[903_940inv;851_885inv]')
---------------------------------------------------------------------------
ParseError                                Traceback (most recent call last)
~/playground/python/hgvs/venv/lib/python3.9/site-packages/hgvs/parser.py in rule_fxn(s)
    129                 try:
--> 130                     return self._grammar(s).__getattr__(rule_name)()
    131                 except ometa.runtime.ParseError as exc:

~/playground/python/hgvs/venv/lib/python3.9/site-packages/parsley.py in invokeRule(*args, **kwargs)
     97                                      [["message", "expected EOF"]], err.trail)
---> 98             raise err
     99         return invokeRule

~/playground/python/hgvs/venv/lib/python3.9/site-packages/ometa/runtime.py in _or(self, fns)
    597                 m = self.input
--> 598                 ret, err = f()
    599                 errors.append(err)

/pymeta_generated_code/pymeta_grammar__Grammar.py in _G_or_212()
    803                 self._trace('| dna_in', (5502, 5510), self.input.position)
--> 804                 _G_apply_213, lastError = self._apply(self.rule_dna_ins, "dna_ins", [])
    805                 self.considerError(lastError, None)

~/playground/python/hgvs/venv/lib/python3.9/site-packages/ometa/runtime.py in _apply(self, rule, ruleName, args)
    494                 memoRec = self.input.setMemo(ruleName,
--> 495                                          [rule(), self.input])
    496             except ParseError as e:

/pymeta_generated_code/pymeta_grammar__Grammar.py in rule_dna_ins(self)
    980                 return (_G_many1_275, self.currentError)
--> 981             _G_consumedby_276, lastError = self.consumedby(_G_consumedby_272)
    982             self.considerError(lastError, 'dna_ins')

~/playground/python/hgvs/venv/lib/python3.9/site-packages/ometa/runtime.py in consumedby(self, expr)
    675         oldInput = self.input
--> 676         _, e = expr()
    677         slice = oldInput.data[oldInput.position:self.input.position]

/pymeta_generated_code/pymeta_grammar__Grammar.py in _G_consumedby_272()
    977                     return (_G_apply_274, self.currentError)
--> 978                 _G_many1_275, lastError = self.many(_G_many1_273, _G_many1_273())
    979                 self.considerError(lastError, None)

/pymeta_generated_code/pymeta_grammar__Grammar.py in _G_many1_273()
    974                     self._trace('<dn', (5909, 5912), self.input.position)
--> 975                     _G_apply_274, lastError = self._apply(self.rule_dna, "dna", [])
    976                     self.considerError(lastError, None)

~/playground/python/hgvs/venv/lib/python3.9/site-packages/ometa/runtime.py in _apply(self, rule, ruleName, args)
    494                 memoRec = self.input.setMemo(ruleName,
--> 495                                          [rule(), self.input])
    496             except ParseError as e:

/pymeta_generated_code/pymeta_grammar__Grammar.py in rule_dna(self)
   2667             self._trace('= dna_iupa', (11894, 11904), self.input.position)
-> 2668             _G_apply_809, lastError = self._apply(self.rule_dna_iupac, "dna_iupac", [])
   2669             self.considerError(lastError, 'dna')

~/playground/python/hgvs/venv/lib/python3.9/site-packages/ometa/runtime.py in _apply(self, rule, ruleName, args)
    494                 memoRec = self.input.setMemo(ruleName,
--> 495                                          [rule(), self.input])
    496             except ParseError as e:

/pymeta_generated_code/pymeta_grammar__Grammar.py in rule_dna_iupac(self)
   2682                 return (_G_python_813, self.currentError)
-> 2683             _G_pred_814, lastError = self.pred(_G_pred_811)
   2684             self.considerError(lastError, 'dna_iupac')

~/playground/python/hgvs/venv/lib/python3.9/site-packages/ometa/runtime.py in pred(self, expr)
    648         if not val:
--> 649             raise e
    650         else:

ParseError:
NM_004006.2:c.940_941ins[903_940inv;851_885inv]
                        ^
Parse error at line 1, column 24: Syntax error. trail: [dna_iupac dna dna_ins dna_edit c_posedit c_variant hgvs_variant]

During handling of the above exception, another exception occurred:

HGVSParseError                            Traceback (most recent call last)
<ipython-input-6-9d483e266162> in <module>
----> 1 var_g = hp.parse_hgvs_variant('NM_004006.2:c.940_941ins[903_940inv;851_885inv]')

~/playground/python/hgvs/venv/lib/python3.9/site-packages/hgvs/parser.py in rule_fxn(s)
    130                     return self._grammar(s).__getattr__(rule_name)()
    131                 except ometa.runtime.ParseError as exc:
--> 132                     raise HGVSParseError("{s}: char {exc.position}: {reason}".format(
    133                         s=s, exc=exc, reason=exc.formatReason()))
    134

HGVSParseError: NM_004006.2:c.940_941ins[903_940inv;851_885inv]: char 24: Syntax error


text:
Converting Control and Case vcf files into a format that can be processed by plink 1.9.0
I have two vcf files, one containing samples with TP53 mutations, and one containing samples with no TP53 mutations. The vcf without the mutations is the control and the vcf with the mutations is the case. How can I combine these two files into something that can be processed by the --assoc function in plink 1.9.0? Do I need to make them into one vcf? Or is it not necessary, and I can just use --assoc on each file after converting each individually?

Both Plink 1.9 [uses --assoc] and Plink 2 [uses --glm] can process VCF files directly; no conversion needed.
I think that Plink will combine multiple files specified on the command line into one file, but if this is not the case then bcftools can be used to merge them.
As described in those input specifications, you can provide Plink with a .fam file (via --fam) to add in phenotype data.
Information about the .fam formats can be found here. It's a 6-column headerless space-separated format containing the following fields:

Family ID ('FID')
Within-family ID ('IID'; cannot be '0')
Within-family ID of father ('0' if father isn't in dataset)
Within-family ID of mother ('0' if mother isn't in dataset)
Sex code ('1' = male, '2' = female, '0' = unknown)
Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)

The phenotype column can be quantitative as well; which is triggered when a numeric value is seen that is not one of the above qualitative codes.

Alright, so should I be able to indicate with a .fam file which samples are from the control group and which ones are from the case group?


text:
Comparing homozygosity of k-mer plots
Attached are two kmer plots from two closely related species. Is that safe to say that the one on the left has higher homozygosity than the one in the right k-mer plot, due to a low to almost flat blue peak (2x) in the left plot? My concern is that the blue peaks in the two plots are not in the same position: left plot at 4n, while right plot at 2n.
I am aware of the grey area in the right plot indicating missing content in assembly, likely due to reads from multiple individuals (correct me if I am wrong).

Model for kmer plot on the left
Formula: y_transform ~ x^transform_exp * length * predict2_0(r1, k, d,
    kmercov, bias, x)

Parameters:
         Estimate Std. Error t value Pr(>|t|)
d       4.019e-02  1.262e-03   31.84   <2e-16 ***
r1      1.014e-02  1.664e-04   60.94   <2e-16 ***
kmercov 1.996e+01  2.452e-02  814.12   <2e-16 ***
bias    1.900e+00  2.495e-02   76.17   <2e-16 ***
length  4.553e+08  8.993e+05  506.34   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4024000 on 1995 degrees of freedom

Number of iterations to convergence: 8
Achieved convergence tolerance: 1.49e-08

Model for kmer plot on the right
Formula: y_transform ~ x^transform_exp * length * predict2_0(r1, k, d,
    kmercov, bias, x)

Parameters:
         Estimate Std. Error t value Pr(>|t|)
d       4.399e-02  1.407e-03   31.27   <2e-16 ***
r1      1.519e-02  1.002e-04  151.64   <2e-16 ***
kmercov 2.069e+01  2.416e-02  856.53   <2e-16 ***
bias    1.697e+00  2.099e-02   80.88   <2e-16 ***
length  5.069e+08  1.056e+06  480.10   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4878000 on 1995 degrees of freedom

Number of iterations to convergence: 7
Achieved convergence tolerance: 1.49e-08

Both the GenomeScope estimated genome sizes are smaller than the assembly (both have microorganisms filtered). The genome size of the left species is ~100Mb smaller than the one on the right, consistently in GenomeScope and the assembly.
GenomeScope Plots:

I would say yes, but it is not that easy to see it on a histogram decomposed to overlapping distributions. While the decomposition to different ditributions is assembly dependent, the overall k-mer histogram is actually directly derived from the sequencing data.
And these k-mer histograms can be used to estimate heterozygosity. For example, using genomescope. You can install the tool locally, or fit the model on a webserver (which seems to be down at the moment, but it should be back on at some point). For more details check the paper. If you get the tool running, you can estimate heterozygosity in the two samples to see how much more/less heterozygous they are.
P.S. I am always thinking about heterozygosity, which is supposed an equivalent of 1 - homozygosity, so I hope the answer will make sense to you.
--- edit ---
so, looking at the gneomescope models and knowing that genome size estimate matches the assembly size, I think it's fair to conclude the right genome is more heterozygous than the left one. However, keep in mind that it is not super great sequencing dataset, as the error peak and heterozygous peak are not clearly separated, hence the estimates of genomic features might be less accurate.


text:
How to find an unidentified splicing variant of a protein?
I have some data indicating there might be a splicing variant of the Arabidopsis Thaliana protein I'm studying that has not been identified. Is there a database of i.e. RNA sequences (transcriptome?) that I could search for it, knowing a part of its sequence?

Please consider accepting the answer below (green tick box). It is comprehensive and help our site stats

There's not going to be a specific database that's particularly thorough for this. The best you can do is one of the following:

Search through TAIR (https://www.arabidopsis.org/servlets/Search?action=new_search&type=protein)
Blast SRA with the full length sequence ( https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&BLAST_PROGRAMS=megaBlast&PAGE_TYPE=BlastSearch&BLAST_SPEC=SRA&SHOW_DEFAULTS=on). This will unfortunately pull up a LOT of related hits for every other related isoform.
Download some public RNA-seq data from SRA/ENA and see if it has supporting data.

You're probably best off doing all of these, where you use blast to find relevant datasets that you then align. My guess is that TAIR won't have everything, since this has probably been seen before but not reported since the experiment wasn't intending to find new proteins.


text:
Removing adapters + primers from fastq files
I am currently in the process of removing adapters and primers from some 16S data acquired, and have a conceptual question regarding this pre-processing step:
So there are a series of different primer sequences (4 for the forward and 4 for the reverse strands), as well as two adapter sequences (1 for each strand) overhanging these primers.
To test this process I am currently using cutadapt. Just to give you an example of what one of my primers vs. adapters look like:
"ACACTCTTTCCCTACACGACGCTCTTCCGATCTGTGCCAGCMGCCGCGGTAA" -> fp
"ACACTCTTTCCCTACACGACGCTCTTCCGATCT" -> a1
So my question is: Given that the adapters are overhanging these primers (a1 is a an exact match with first N base pairs in fp), could we just use cutadapt once, in order to remove the primers + adapter sequences in one go?


text:
Length of Contigs in Transcriptome and Whole Genome Assembly
Why are there shorter contigs from transcriptome assembly than from a whole genome assembly? I know the difference between transcriptome and genome, but don't really understand what contigs are in the context of sequencing in bioinformatics, and why are the contigs shorter in transcriptome assembly than in whole genome assembly. Anyone mind explaining? Much appreciated.

The longest possible transcriptome contig reflects the longest possible transcript. In the human genome, that's possibly Titin, which is ~35K aa (that's ~105Kb) and will be longer pre-splicing.
The longest possible genome contig reflects the longest chromosome. In the human genome that's chr1 at ~250Mb.
In practise your contigs/scaffolds will be a lot shorter for various reasons (e.g. hard to sequence the telomeres and repeats) but it illustrates the difference since all the chromosomes consist of contiguous DNA but many separate RNA transcripts (and not all DNA is transcribed).


text:
Visualization of mutations in (part of) viral genome - existing tools?
I'm trying to visualize mutations in SARS-CoV-2 sequences, somewhat like it's done on outbreak.info:

but for individual sequences, probably on the base of the vcf files generated by the ARTIC pipeline. A pipeline focused on virus nanopore output.
The ready-made tools I've found for visualization of mutations so far were specifically tailored to human sequences and/or for visualizing aggregated data for multiple genomes, or were web tools, which can't be used for medical data due to confidentiality.
Is there any existing offline tool that

is Linux-based
can be integrated into an automated pipeline
is tailored towards viral genomes or is "organism-agnostic"
is based on comparing an individual sequence to a reference sequence
shows the sequence as "blocks" of genes with the mutations highlighted
(ideally plays well with MultiQC for integration into an existing report, but that's a nice-to-have)

?

I ended up using lolliplots from the trackViewer R package; it can do VCF-based visualization and is flexible enough I can feed it data from other sources (I wound up parsing Nextclade mutation calls for the amino acids since I had these anyway).

Thanks for your update, this really useful. The Vignette package I do recognise, it does look really cool


text:
List of all Latin word parts used in Biological Taxonomy?
I seem to have determined that there are 287,614 unique Latin taxonomic names used in the biological taxonomy according to the Integrated Taxonomic Information System (ITIS) taxonomy taxon.txt file. First of all, that is an enormous number of unique names, I would have hoped they could be reused to a greater degree but I guess not. There are 90k some odd genera names, the rest are species names, or certain names might be duplicated throughout the complete taxonomy.
What I would like to do is take the taxonomy Latin names and convert them into their English equivalents. Not their "common names", as that can't be done by looking at the Latin name alone. But by taking their Latin name and making it into the equivalent English literal meaning. Wikipedia has a nice list on commonly used taxonomic affixes, and another on Latin and Greek words commonly used in systematic names, however a quick glance at the list and they both are missing -bacter- as one infix, for a single example, or vorans, or many many others. The Wikipedia list seems to only contain a small fraction of the names used in the actual complete taxonomy.
What I'm wondering is if there exist anywhere online in text format (so it can be copy-pasted) a complete list of all the Latin affixes (prefixes, infixes, suffixes, etc.) so I can accomplish my goal. The goal being, programmatically convert the Latin names to their literal English equivalent. I would split apart the Latin name based on some database of Latin affixes, and then plop out the English word matching the affix. It should take only a few seconds for a computer to do this transformation, assuming the list mapping Latin names to English names exists already.
Does anything like this exist anywhere on the web?
This calls itself the "worlds largest affix directory", yet it still doesn't have these common ones... Here's another Wikipedia page, perhaps this is as good as it is going to get?

The keyword is calque: https://en.wikipedia.org/wiki/Calque
This question is very much a question of automated calque translation and may be more suited for a linguistics or vanilla Stack Overflow (as the bio part isn't relevant). A resource that may be doing what you are after already is LPSN at https://www.bacterio.net/: This has the etymology of properly described bacteria taxa, including humorous notes about bad grammar... Hello Pelagibacter ubique.

I anticipate two problems with this project - many of the latin names are just names of biologists, therefore non-meaningful words in latin. Second, latin names are usually give by people that don't really know latin, therefore there are quite often grammatical errors (e.g. Timema monikensis/ense)

Maybe you can use NCBI Taxonomy instead, which is a full hierarchical taxonomy with >700K taxa. You can use Entrez utilities to programmatically query NCBI taxonomy to retrieve Latin or common names as well as other metadata.

This doesn't answer the question :/ How specifically do I get the latin name literal English meaning from it though?


text:
Complete List of IUPAC Chemical Naming Conventions for Formal and Trivial Names?
Two-parted question:

Where is a cheatsheet or complete list of the "formal" IUPAC naming conventions?
Where is a cheatsheet or complete list of the "trivial" or "informal" chemical names?

Does the Wikipedia page have all of the naming conventions, or are some missing? Then this page titled "Trivial and semisystematic names retained for naming organic compounds" seems to link to some of the trivial names, but it is in a terrible format for a computer to read (I would like a spreadsheet or some sort of structured data format so I can use these in code) and I'm not sure it's complete. I have searched around on the web but don't seem to find anything that is complete and ideally as structured data. Does anything like this exist? If not structured data, what is the easiest document to read that is complete that I could manually convert into structured data? Also, any of this should ideally be a free resource.

IUPAC Reports
The International Union of Pure and Applied Chemistry (IUPAC) produces recommendations on nomenclature for chemical compounds. They produce complete reports, and brief guides. However, creating a systematic name from an arbitrary structure is not simple.
It probably depends what 'space' of structures you want to deal with. For simple organic compounds (branched, fused rings, non-bridged) it might be feasible, but there are many edge cases and tricky topics, such as stereochemistry, ligand priority ordering, multiple bridges, and so on.
For example, this is the 'Blue Book' : https://iupac.qmul.ac.uk/BlueBook/ which is a set of PDF docs. Converting these to formal rules would be time consuming, to say the least!
Informal Names
A resource like ChemSpider provides multiple informal names for compounds, with an API to retrieve data. Seems like (from a quick test) you can get the 'CommonName'  field from the details. For example, for a random recordId:
https://api.rsc.org/compounds/v1/records/325/details
you get this return value:
{"id": 325,"commonName": "2-Methoxy-4-vinylphenol"}
so it gives you one common name. Seemingly also for benzene (id 236).


text:
Assigning case and control label to different samples of large merged file to input into a plink 1.9.0 association test
I have a large vcf file with both case and control samples in the file. I am planning to input the vcf into the --assoc function of plink with the --fam parameter that contains a .fam file that specifies which samples are case and which are control. I want to make a .fam file with the case and control samples labelled in the phenotype value column. How can I assign case/control to different samples for the .fam of my merged vcf? My big vcf file doesn't indicate whether the samples are case or control. I have two tsv files, one with a list of control samples and one with a list of case samples. Can I use these two files to specify in the .fam file which samples are case and which are control?
Could I also set the specific gene mutation as the case and samples without the mutation as the control, and run an association test based on those parameters?

It would help to see what your control/case tsv files look like. Here is the .fam file specification, IDK if you have family information based on what you've posted so many of the entries may simply be zero.

See this post for some information about preparing files: https://bioinformatics.stackexchange.com/questions/15385/making-plink-compatible-files-from-vcf-file-without-phenotype-information. If you have comfort with scripting tools such as python, you could just update the .fam file from your TSV files given that phenotype info is missing from your VCF.


text:
MarkDuplicatesSpark failing with cryptic error message. MarkDuplicates succeeds
I have been trying to follow the GATK Best Practice Workflow for 'Data pre-processing for variant discovery' (https://gatk.broadinstitute.org/hc/en-us/articles/360035535912).
This has all been run on Windows Subsystem for Linux 2 on the Bash shell.
I started off with FASTQ files from IGSR (https://www.internationalgenome.org/data-portal) and performed alignment with Bowtie2 (instead of BWA).
This produced a SAM file which I converted to a BAM file using SAMTOOLS.
I have been trying to use MarkDuplicatesSpark (https://gatk.broadinstitute.org/hc/en-us/articles/4409897162139-MarkDuplicatesSpark) but am getting an output that I can't make sense of.
I am running the following command:
gatk MarkDuplicatesSpark -I HG00102_hGRCH38_exome_aignment.bam -O HG00102_hGRCH38_exome_aignment.marked_duplicates.bam --remove-sequencing-duplicates --spark-master local[*]

I am getting the following output:
Using GATK jar /mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar MarkDuplicatesSpark -I HG00102_hGRCH38_exome_aignment.bam -O HG00102_hGRCH38_exome_aignment.marked_duplicates.bam --remove-sequencing-duplicates --spark-master local[*]

22:45:08.487 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
Dec 19, 2021 10:45:08 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
22:45:08.748 INFO  MarkDuplicatesSpark - ------------------------------------------------------------
22:45:08.748 INFO  MarkDuplicatesSpark - The Genome Analysis Toolkit (GATK) v4.2.4.0
22:45:08.748 INFO  MarkDuplicatesSpark - For support and documentation go to https://software.broadinstitute.org/gatk/
22:45:08.748 INFO  MarkDuplicatesSpark - Executing as *removed* on Linux v4.19.128-microsoft-standard amd64
22:45:08.749 INFO  MarkDuplicatesSpark - Java runtime: OpenJDK 64-Bit Server VM v11.0.9.1+1-Ubuntu-0ubuntu1.18.04
22:45:08.749 INFO  MarkDuplicatesSpark - Start Date/Time: December 19, 2021 at 10:45:08 PM GMT
22:45:08.749 INFO  MarkDuplicatesSpark - ------------------------------------------------------------
22:45:08.749 INFO  MarkDuplicatesSpark - ------------------------------------------------------------
22:45:08.750 INFO  MarkDuplicatesSpark - HTSJDK Version: 2.24.1
22:45:08.750 INFO  MarkDuplicatesSpark - Picard Version: 2.25.4
22:45:08.751 INFO  MarkDuplicatesSpark - Built for Spark Version: 2.4.5
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.COMPRESSION_LEVEL : 2
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
22:45:08.751 INFO  MarkDuplicatesSpark - Deflater: IntelDeflater
22:45:08.751 INFO  MarkDuplicatesSpark - Inflater: IntelInflater
22:45:08.751 INFO  MarkDuplicatesSpark - GCS max retries/reopens: 20
22:45:08.752 INFO  MarkDuplicatesSpark - Requester pays: disabled
22:45:08.752 INFO  MarkDuplicatesSpark - Initializing engine
22:45:08.752 INFO  MarkDuplicatesSpark - Done initializing engine
Using Spark's default log4j profile: org/apache/spark/log4j-defaults.properties
21/12/19 22:45:09 WARN Utils: Your hostname, *hostname removed* resolves to a loopback address: *IP removed*; using 172.27.147.65 instead (on interface eth0)
21/12/19 22:45:09 WARN Utils: Set SPARK_LOCAL_IP if you need to bind to another address
WARNING: An illegal reflective access operation has occurred
WARNING: Illegal reflective access by org.apache.spark.unsafe.Platform (file:/mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar) to method java.nio.Bits.unaligned()
WARNING: Please consider reporting this to the maintainers of org.apache.spark.unsafe.Platform
WARNING: Use --illegal-access=warn to enable warnings of further illegal reflective access operations
WARNING: All illegal access operations will be denied in a future release
21/12/19 22:45:10 INFO SparkContext: Running Spark version 2.4.5
21/12/19 22:45:10 WARN NativeCodeLoader: Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
21/12/19 22:45:10 INFO SparkContext: Submitted application: MarkDuplicatesSpark
21/12/19 22:45:10 INFO SecurityManager: Changing view acls to: abgane
21/12/19 22:45:10 INFO SecurityManager: Changing modify acls to: abgane
21/12/19 22:45:10 INFO SecurityManager: Changing view acls groups to:
21/12/19 22:45:10 INFO SecurityManager: Changing modify acls groups to:
21/12/19 22:45:10 INFO SecurityManager: SecurityManager: authentication disabled; ui acls disabled; users  with view permissions: Set(abgane); groups with view permissions: Set(); users  with modify permissions: Set(abgane); groups with modify permissions: Set()
21/12/19 22:45:10 INFO Utils: Successfully started service 'sparkDriver' on port 35319.
21/12/19 22:45:10 INFO SparkEnv: Registering MapOutputTracker
21/12/19 22:45:11 INFO SparkEnv: Registering BlockManagerMaster
21/12/19 22:45:11 INFO BlockManagerMasterEndpoint: Using org.apache.spark.storage.DefaultTopologyMapper for getting topology information
21/12/19 22:45:11 INFO BlockManagerMasterEndpoint: BlockManagerMasterEndpoint up
21/12/19 22:45:11 INFO DiskBlockManager: Created local directory at /tmp/blockmgr-1e5e323f-dd90-4093-980c-0ceda43349c2
21/12/19 22:45:11 INFO MemoryStore: MemoryStore started with capacity 3.6 GB
21/12/19 22:45:11 INFO SparkEnv: Registering OutputCommitCoordinator
21/12/19 22:45:11 INFO Utils: Successfully started service 'SparkUI' on port 4040.
21/12/19 22:45:11 INFO SparkUI: Bound SparkUI to 0.0.0.0, and started at http://172.27.147.65:4040
21/12/19 22:45:11 INFO Executor: Starting executor ID driver on host localhost
21/12/19 22:45:11 INFO Utils: Successfully started service 'org.apache.spark.network.netty.NettyBlockTransferService' on port 33189.
21/12/19 22:45:11 INFO NettyBlockTransferService: Server created on 172.27.147.65:33189
21/12/19 22:45:11 INFO BlockManager: Using org.apache.spark.storage.RandomBlockReplicationPolicy for block replication policy
21/12/19 22:45:11 INFO BlockManagerMaster: Registering BlockManager BlockManagerId(driver, 172.27.147.65, 33189, None)
21/12/19 22:45:11 INFO BlockManagerMasterEndpoint: Registering block manager 172.27.147.65:33189 with 3.6 GB RAM, BlockManagerId(driver, 172.27.147.65, 33189, None)
21/12/19 22:45:11 INFO BlockManagerMaster: Registered BlockManager BlockManagerId(driver, 172.27.147.65, 33189, None)
21/12/19 22:45:11 INFO BlockManager: Initialized BlockManager: BlockManagerId(driver, 172.27.147.65, 33189, None)
22:45:11.804 INFO  MarkDuplicatesSpark - Spark verbosity set to INFO (see --spark-verbosity argument)
21/12/19 22:45:11 INFO GoogleHadoopFileSystemBase: GHFS version: 1.9.4-hadoop3
21/12/19 22:45:12 INFO MemoryStore: Block broadcast_0 stored as values in memory (estimated size 172.9 KB, free 3.6 GB)
21/12/19 22:45:13 INFO MemoryStore: Block broadcast_0_piece0 stored as bytes in memory (estimated size 35.5 KB, free 3.6 GB)
21/12/19 22:45:13 INFO BlockManagerInfo: Added broadcast_0_piece0 in memory on 172.27.147.65:33189 (size: 35.5 KB, free: 3.6 GB)
21/12/19 22:45:13 INFO SparkContext: Created broadcast 0 from newAPIHadoopFile at PathSplitSource.java:96
21/12/19 22:45:13 INFO MemoryStore: Block broadcast_1 stored as values in memory (estimated size 172.9 KB, free 3.6 GB)
21/12/19 22:45:13 INFO MemoryStore: Block broadcast_1_piece0 stored as bytes in memory (estimated size 35.5 KB, free 3.6 GB)
21/12/19 22:45:13 INFO BlockManagerInfo: Added broadcast_1_piece0 in memory on 172.27.147.65:33189 (size: 35.5 KB, free: 3.6 GB)
21/12/19 22:45:13 INFO SparkContext: Created broadcast 1 from newAPIHadoopFile at PathSplitSource.java:96
21/12/19 22:45:13 INFO FileInputFormat: Total input files to process : 1
21/12/19 22:45:13 INFO SparkUI: Stopped Spark web UI at http://172.27.147.65:4040
21/12/19 22:45:13 INFO MapOutputTrackerMasterEndpoint: MapOutputTrackerMasterEndpoint stopped!
21/12/19 22:45:13 INFO MemoryStore: MemoryStore cleared
21/12/19 22:45:13 INFO BlockManager: BlockManager stopped
21/12/19 22:45:13 INFO BlockManagerMaster: BlockManagerMaster stopped
21/12/19 22:45:13 INFO OutputCommitCoordinator$OutputCommitCoordinatorEndpoint: OutputCommitCoordinator stopped!
21/12/19 22:45:13 INFO SparkContext: Successfully stopped SparkContext
22:45:13.851 INFO  MarkDuplicatesSpark - Shutting down engine
[December 19, 2021 at 10:45:13 PM GMT] org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark done. Elapsed time: 0.09 minutes.
Runtime.totalMemory()=419430400
java.lang.IllegalArgumentException: Unsupported class file major version 55
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:166)
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:148)
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:136)
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:237)
    at org.apache.spark.util.ClosureCleaner$.getClassReader(ClosureCleaner.scala:49)
        at org.apache.spark.util.FieldAccessFinder$$anon$3$$anonfun$visitMethodInsn$2.apply(ClosureCleaner.scala:517)
        at org.apache.spark.util.FieldAccessFinder$$anon$3$$anonfun$visitMethodInsn$2.apply(ClosureCleaner.scala:500)
        at scala.collection.TraversableLike$WithFilter$$anonfun$foreach$1.apply(TraversableLike.scala:733)
    at scala.collection.mutable.HashMap$$anon$1$$anonfun$foreach$2.apply(HashMap.scala:134)
    at scala.collection.mutable.HashMap$$anon$1$$anonfun$foreach$2.apply(HashMap.scala:134)
    at scala.collection.mutable.HashTable$class.foreachEntry(HashTable.scala:236)
        at scala.collection.mutable.HashMap.foreachEntry(HashMap.scala:40)
        at scala.collection.mutable.HashMap$$anon$1.foreach(HashMap.scala:134)
    at scala.collection.TraversableLike$WithFilter.foreach(TraversableLike.scala:732)
    at org.apache.spark.util.FieldAccessFinder$$anon$3.visitMethodInsn(ClosureCleaner.scala:500)
    at org.apache.xbean.asm6.ClassReader.readCode(ClassReader.java:2175)
    at org.apache.xbean.asm6.ClassReader.readMethod(ClassReader.java:1238)
    at org.apache.xbean.asm6.ClassReader.accept(ClassReader.java:631)
    at org.apache.xbean.asm6.ClassReader.accept(ClassReader.java:355)
    at org.apache.spark.util.ClosureCleaner$$anonfun$org$apache$spark$util$ClosureCleaner$$clean$14.apply(ClosureCleaner.scala:307)
    at org.apache.spark.util.ClosureCleaner$$anonfun$org$apache$spark$util$ClosureCleaner$$clean$14.apply(ClosureCleaner.scala:306)
        at scala.collection.immutable.List.foreach(List.scala:392)
        at org.apache.spark.util.ClosureCleaner$.org$apache$spark$util$ClosureCleaner$$clean(ClosureCleaner.scala:306)
    at org.apache.spark.util.ClosureCleaner$.clean(ClosureCleaner.scala:162)
        at org.apache.spark.SparkContext.clean(SparkContext.scala:2326)
        at org.apache.spark.SparkContext.runJob(SparkContext.scala:2100)
        at org.apache.spark.SparkContext.runJob(SparkContext.scala:2126)
        at org.apache.spark.rdd.RDD$$anonfun$collect$1.apply(RDD.scala:990)
    at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:151)
    at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:112)
    at org.apache.spark.rdd.RDD.withScope(RDD.scala:385)
    at org.apache.spark.rdd.RDD.collect(RDD.scala:989)
    at org.apache.spark.RangePartitioner$.sketch(Partitioner.scala:309)
    at org.apache.spark.RangePartitioner.<init>(Partitioner.scala:171)
    at org.apache.spark.RangePartitioner.<init>(Partitioner.scala:151)
    at org.apache.spark.rdd.OrderedRDDFunctions$$anonfun$sortByKey$1.apply(OrderedRDDFunctions.scala:62)
        at org.apache.spark.rdd.OrderedRDDFunctions$$anonfun$sortByKey$1.apply(OrderedRDDFunctions.scala:61)
        at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:151)
        at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:112)
        at org.apache.spark.rdd.RDD.withScope(RDD.scala:385)
        at org.apache.spark.rdd.OrderedRDDFunctions.sortByKey(OrderedRDDFunctions.scala:61)
        at org.apache.spark.api.java.JavaPairRDD.sortByKey(JavaPairRDD.scala:936)
        at org.broadinstitute.hellbender.utils.spark.SparkUtils.sortUsingElementsAsKeys(SparkUtils.java:165)
        at org.broadinstitute.hellbender.utils.spark.SparkUtils.sortReadsAccordingToHeader(SparkUtils.java:143)
        at org.broadinstitute.hellbender.utils.spark.SparkUtils.querynameSortReadsIfNecessary(SparkUtils.java:306)
        at org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark.mark(MarkDuplicatesSpark.java:206)
        at org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark.mark(MarkDuplicatesSpark.java:270)
        at org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark.runTool(MarkDuplicatesSpark.java:350)
        at org.broadinstitute.hellbender.engine.spark.GATKSparkTool.runPipeline(GATKSparkTool.java:546)
        at org.broadinstitute.hellbender.engine.spark.SparkCommandLineProgram.doWork(SparkCommandLineProgram.java:31)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:140)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
        at org.broadinstitute.hellbender.Main.main(Main.java:289)
21/12/19 22:45:13 INFO ShutdownHookManager: Shutdown hook called
21/12/19 22:45:13 INFO ShutdownHookManager: Deleting directory /tmp/spark-f49df209-c946-40e3-96ca-9d8125c3a733

As an alternative, I am trying to sort the BAM file and use MarkDuplicates instead (https://gatk.broadinstitute.org/hc/en-us/articles/4409924785691-MarkDuplicates-Picard-).This appears to work fine. I obviously want to make use of multiprocessing though so ideally would like to get the Spark implementation going.
Does anyone know what might be causing this command to fail?
Any ideas for a fix?
Any alternative software that might be useful?
Thanks for your time,
Angus

I also got a spark error when running MarkDuplicatesSpark. however, when I look into the error message, I found mine is due to network time out: 21/12/21 18:08:34 WARN HeartbeatReceiver: Removing executor driver with no recent heartbeats: 134158 ms exceeds timeout 120000 ms 21/12/21 18:08:37 ERROR TaskSchedulerImpl: Lost executor driver on localhost: Executor heartbeat timed out after 134158 ms.  I was able to fix it by adding --conf 'spark.network.timeout=1000000' to the command line. it looks like yours is probably due to another issue.

You need to set your Java version to Java 8.
If you are using conda, run this in your conda environment:
conda install openjdk=8


text:
What are Dummy fastqs?
I had to do a file transfer for in our system to run from fastqs and I came across the term "dummy fastq", not sure what exactly it means and what is the purpose of it in the workflow. Can someone answer?

Probably some made-up data solely for the sake of testing a workflow but actually impossible to know what random files are sitting on your servers.


text:
Using different "principal components" in UMAP visualizations
In certain UMAP visualization, it appears that the visualizations show the first UMAP "principal component" plotted against the second. For example, like this:

Apparently, this lets us use the first and second UMAP "principal components" to tell apart the different cell clusterings in a dataset. I'm curious about if we were to use e.g. the second and the third UMAP "principal components" instead - would this let us tell apart different cell clusterings according to the second and third features that explain those clusterings' largest sources of variation? If someone could explain this to me, I would be very grateful, and if you would like to correct whatever misunderstandings I have about UMAP as well, that would be great too - I'd be the first one to tell you that I don't understand it well!

The clustering is just visualized onto the UMAP, it is based on the principal components of a PCA, don't mix that up. UMAP in many cases is just a 2D compression of the PCA results for visualization. Nice to have and good to get an idea how data look, but that's all about it.

You can plot UMAP3 too in 3 dimensions, but that's about the limit for viz. You could plot UMAP2 vs. UMAP3, UMAP3 vs. UMAP4 etc. but I'm not sure how informative it is. Instead, you can focus on other analyses like clustering, differential expression, cell type annotation, marker gene analysis, and other machine learning approaches etc. Ultimately, UMAP/t-SNE are nice tools for viz but they can be over interpreted- it is inherently impossible to summarise very highly dimensional data in just a few dimensions without losing information- see here for an (overly) provocative take.


text:
Remove variable sequence component within a tree text file
I have a gene tree file of 436 orthologue genes from 6 species. I want to remove unwanted extensions as it looks massy after visualization. My file looks like:
(TRINITY_Clupea_DN5452_c0_g1_i1.p1:0.0824467436,TRINITY_Engraulis_DN43599_c0_g1_i1.p1:0.1634781085)100:0.0876433106,TRINITY_Sardina_DN15766_c0_g1_i2.p1:0.0164132018)..................
What i need:
(Clupea_DN5452:0.0824467436,Engraulis_DN43599:0.1634781085)100:0.0876433106,Sardina_DN15766:0.0164132018)..................
As "TRINITY" is identicial, i can remove it using sed. But after the species name ids are not identical. And i only need the 2nd and 3rd part of the identifier.
It will be helpful for me if i get some suggestion. Thanks.

Lots of ways to do this, but here's one way using sed:
sed 's/TRINITY_\([^_]\+\)_\([^_]\+\)_[^:]\+/\1_\2/g' file.tree

The above applies a substitution to capture the first two fields using parentheses (which need to be escaped). Each field captures any character not an underscore one or more times (using the "+" symbol, which also needs to be escaped). Finally, we also match any character not a colon one or more times. The replacement string can be just the first two capture groups separated with an underscore.

Hi Steve,
Thanks for your warm response. But my file is a ".tree" file, not a fasta file since the file dosent start with a ">" symbol. Like as following the file looks:
(TRINITY_Clupea_DN5452_c0_g1_i1.p1:0.0824467436,TRINITY_Engraulis_DN43599_c0_g1_i1.p1:0.1634781085)100:0.0876433106,TRINITY_Sardina_DN15766_c0_g1_i2.p1:0.0164132018)..................
And i want like:
(Clupea_DN5452:0.0824467436,Engraulis_DN43599:0.1634781085)100:0.0876433106,Sardina_DN15766:0.0164132018)..................
Want to keep only the 2nd and 3rd part from the id.
Thanks.

If you want to remove component/gene/isoform parts from a Trinity ID, you can do s/_c[0-9]+_.*?:/:/ to remove anything between the first c with only numbers after it up to the next colon. Combining this with the TRINITY removal at the start, it'd be something like this:
sed 's/TRINITY_\(.*?\)_c[0-9]+_.*?:/\1:/g' file.tree

I have tried like as follows:
sed 's/TRINITY_(.?)c[0-9]+.?:/\1:/g' OG0000004.tre
But is priniting like as:
(TRINITY_Engraulis_DN849_c0_g1_i1.p1:0.05941030350315353731,(((((TRINITY_Tenualosa_DN4978_c0_g1_i1.p1:0.07382876155314552213,TRINITY_Alosa_DN5385_c0_g1_i1.p1:0.01287649339098437609):0.00905926786283233615,TRINITY_Sardina_DN264_c0_g1_i7.p1:0.02335748047027334101):0.04151957750042657019,((TRINITY_Colia_DN12738_c0_g1_i1.p1:0.04097580399145889923,(TRINITY_Engraulis_DN8273_c0_g2_i3.p1:0.04815781217046399831,(((((TRINITY_Engraulis_DN16258_c3_g1_i1.p1:0.08818477156986732901,....

@kishor-kumar-sarker Looks like you didn't escape the parentheses


text:
Why does randomized motif search work?
I'm looking for intuition for why a randomized motif search works. My current thinking is as follows:
We are selecting many random kmers from our DNA sequences. The chosen kmers will bias the profile matrix to selecting kmers like them.
Given any particular k-mer chosen, there are two possibilities:

We've selected a meaningless kmer, so the profile will bias selecting random (ie high entropy) kmers and the score will stop improving soon.
We've captured a portion of the true motif. Since the motif isn't random (ie doesn't have high entropy) subsequent iterations using this profile matrix will pick motifs closer to the real motif and the score will keep improving.

Since we give this algorithm many tries and at some point it will eventually arrive at possibility 2 and come to a good conclusion.
Is this a correct understanding of why randomized motif search works? Is there a more rigorous way to prove this?

Update
With regard to comment, the proof of effectiveness of Monte Carlo algorithms like this is a little thorny- if you look at the wikipedia page, you will see that a proof of statistical consistency for any Monte Carlo algorithms took a long time:

From 1950 to 1996, all the publications on Sequential Monte Carlo methodologies, including the pruning and resample Monte Carlo methods introduced in computational physics and molecular chemistry, present natural and heuristic-like algorithms applied to different situations without a single proof of their consistency, nor a discussion on the bias of the estimates and on genealogical and ancestral tree based algorithms. The mathematical foundations and the first rigorous analysis of these particle algorithms were written by Pierre Del Moral in 1996.

Here is a possibly more general reference by Del Moral in 1998. The paper is somewhat dense and I am not an expert in the field so I don't think I'll do it justice, but I think that we could probably figure out a similar mathematical proof for random motif search given time and the right problem formulation.
I sense that your concern might be with Monte Carlo algorithms generally rather than random motif search specifically, so hopefully that gets close to your question.
It does look like people have done heuristic analysis of similar algorithms, but they are not formal proofs.
Original answer
I believe that you are more or less correct.
For example, this presentation's walkthrough of the algorithm (slides 35-36) specifically refers to Greedy randomized profile motif searches. By making the greediness of the algorithm explicit, we can see that there is no claim of finding an optimal solution in any particular iteration, thus the need for multiple retries to get somewhere close to approximate solution.
One additional point made in that presentation is that this algorithm works best when coupled to approaches such as Gibbs sampling for targeting the sampled k-mer chain in regions with higher scores (less entropic k-mers).
A different learning resource with similar content (and a little more intuition) is available here.

That's actually the course I'm taking haha.
I guess what I'm wondering is if we can prove this stuff rigorously. It makes intuitive sense and I'm glad I got the intuition right but I do wonder how we know our intuition is correct.

@Moo see update- I think that addresses the Monte Carlo question, but I can't tell if anyone has done a similar proof for this algo specifically.

Thank you for the update. I think you're right that the root of my question is why Monte Carlo algorithms work.