text:
RagTag patch error--"Tuple index out of range"
This question was also asked on GitHub
I'm trying to correct a long-read assembly with a short-read scaffold; I'm hoping to fill in the short gaps in the scaffold with the matching long-read sections. The RagTag "patch" function, if I understand correctly, is supposed to do exactly this. I use this command:
ragtag.py patch --fill-only liberibacter_central_ragtag.fasta liberibacter_central_assembly.fasta

But when I attempt to run it I get the following error:
Tue Feb  7 09:55:17 2023 --- VERSION: RagTag v2.1.0
Tue Feb  7 09:55:17 2023 --- WARNING: This is a beta version of `ragtag patch`
Tue Feb  7 09:55:17 2023 --- CMD: ragtag.py patch --fill-only liberibacter_central_assembly.fasta liberibacter_central_ragtag.fasta
Tue Feb  7 09:55:17 2023 --- WARNING: Without '-u' invoked, some component/object AGP pairs might share the same ID. Some external programs/databases don't like this. To ensure valid AGP format, use '-u'.
Tue Feb  7 09:55:17 2023 --- INFO: Running: ragtag_splitasm.py -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.agp /xdisk/baltrus/moshesteyn/liberibacter_central_assembly.fasta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Finished running : ragtag_splitasm.py -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.agp /xdisk/baltrus/moshesteyn/liberibacter_central_assembly.fasta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Running: ragtag_rename.py /xdisk/baltrus/moshesteyn/liberibacter_central_ragtag.fasta -p qseq -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.agp > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Finished running : ragtag_rename.py /xdisk/baltrus/moshesteyn/liberibacter_central_ragtag.fasta -p qseq -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.agp > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Writing: /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.comps.fasta
Tue Feb  7 09:55:17 2023 --- INFO: Mapping the query genome to the target genome
Tue Feb  7 09:55:17 2023 --- INFO: Running: nucmer --maxmatch -l 100 -c 500 -p /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta.log
Tue Feb  7 09:55:19 2023 --- INFO: Finished running : nucmer --maxmatch -l 100 -c 500 -p /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta.log
Tue Feb  7 09:55:19 2023 --- INFO: Running: ragtag_delta2paf.py /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.paf 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:19 2023 --- INFO: Finished running : ragtag_delta2paf.py /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.paf 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:19 2023 --- INFO: Reading whole genome alignments
Tue Feb  7 09:55:19 2023 --- INFO: Filtering and merging alignments
Tue Feb  7 09:55:19 2023 --- INFO: Building a scaffold graph from the contig AGP file
Tue Feb  7 09:55:19 2023 --- INFO: Building a scaffold graph from the target/query mappings
Tue Feb  7 09:55:19 2023 --- INFO: Computing a matching solution to the scaffold graph
Traceback (most recent call last):
File "/home/u30/moshesteyn/.conda/envs/py382/bin/ragtag_patch.py", line 518, in <module>
    main()
File "/home/u30/moshesteyn/.conda/envs/py382/bin/ragtag_patch.py", line 497, in main
    match_psg = aln_psg.max_weight_matching()
File "/home/u30/moshesteyn/.conda/envs/py382/lib/python3.8/site-packages/ragtag_utilities/ScaffoldGraph.py", line 824, in max_weight_matching
    for cc in nx.connected_components(G=cover_graph):
File "/home/u30/moshesteyn/.conda/envs/py382/lib/python3.8/site-packages/networkx/classes/backends.py", line 134, in wrapper
    graph = args[0]
IndexError: tuple index out of range

Does anyone know what might be causing this error? No one else seems to have run into it on the github. I've run it with both files individually against a known reference, as well as reversing the order of the files and removing the "fill-only" option, and I have the same problem. I've looked at several of the output files, and they seem normal, at least in that they aren't empty or filled with zeroes.

Were any errors/warnings logged in any of the *.err or *.log files (in the ragtag_output folder)?

The output mentions potential AGP issues if not using -u; have you tried running it with -u included?

Can you post some assembly stats for your various assemblies? I could imagine that a filtering process under the hood removes some small artefactual contigs and exposes some sort of edge case, but it's not totally clear from what you've shown so far. For example, do you have any tiny short-read contigs that are all GGGG or something?

@Steve no errors noted in the .log, and .err is empty.

-u makes no difference, and the shortest contigs are non-repetitive between 190 and 600 bp. What kinds of stats are you looking for--spades log?

Sorry to spam the comments. Attempting to run my own patch, I've mapped the contigs individually to the reference, but for two short ones Bowtie doesn't find matches at all. Might this be causing the problem?

@schmiggle you will need to use @ to tag us when you comment, or we won't see it. Yes, any sort of stats on the assembly would be helpful, spades/quast would be fine. I doubt that the issue is a failure to find a match as that is likely a common case. I'd suggest simply filtering out contigs <500bp because they are unlikely to be very informative.


text:
how to do manual search on NCBI with biopython and add link to publication
I'm trying to reproduce a complex manual search on NCBI using BioPython, but obviously something is not right. My command line is not working, but I'd especially like a way to have an accurate search where I can select only the genomes related to a published paper. Is there an easy way to do that?
I want to limit my results to only those genomes that i) have an NCBI BioProject entry and ii) that BioProject entry contains a Publications field.
For example, let's say I want all the genomes of Microbacterium, the latest version, from RefSeq and not partial (this is what I tried to write but something is not working) and I would like to add one more filter. Take, among the remaining/selected genomes, only the ones related to published papers.
For this strain, we have a related paper:
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA669461/
For this one, we do not:
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA608230/
and I really would like to not have to check 100 genomes manually, but I don't know how to manage it.
I tried:
query = "Microbacterium[Organism] AND latest refseq[filter] AND (all[filter] NOT
partial[filter])"
handle = Entrez.esearch(term=query, db="Assembly", retmax=600)
ids = Entrez.read(handle)["IdList"]

And after this, I'm just using a for loop to get genomes, but my query is not good and it is not handling the problem of the publication.
Thanks for any help!

Edit : as I need to search in the db Assembly unfortunately the answer with the flag author or date is not working

Entrez has the tags [journal], [volume], [page number], [title], [year] and [author] to identify all sequences from a given publication. You will find a full list of tags here from the drop down menu.
Thus,
query = "Microbacterium[Organism] AND Nature[Journal] AND 528[volume] AND Bai[author] AND 364[page number]"

This will identify all Microbacterium genome sequences associated with Bai et al (2015) Nature 528:364- and is here.
For practical purposes I routinely search using the organism and author tags for their sequence output.

The example given (post-edit) was

Kim KK et al., "Microbacterium aquimaris sp. nov., isolated from seawater.", Int J Syst Evol Microbiol, 2008 Jul;58(Pt 7):1616-20

query = "Microbacterium[Organism] AND 58[volume] AND Kim[author] AND 1616[page number] AND 2008[year]"

I removed the [journal] tag because the abbreviation it uses is a bit unusual and replaced it with [year]. Still works here.

From the comments ... if only published sequences were required then using the range function and [publication date] is one way in ...
query = "Microbacterium[Organism] AND 2000/1/1:2020/1/1[publication date]"

The best strategy is to use the [author] because data production tends to be dominated by a small number of groups.

Comments of Comments
Why isn't searching working. The following search works
query = "Microbacterium[Organism] AND (all[filter] NOT partial[filter])"

The term it doesn't like is latest refseq[filter] if that is removed it will work. refseq works, but it doesn't like latest.

I think the OP wants to limit the results to only those genomes whose NCBI BioProject entry includes a Publications field.

@terdon I don't believe there is a tag [publication], there is [publication date], which is better than [year]. However, it's the same answer, the sequence output of most fields are represented only by a small number of authors and only those authors are searched.

@M__ thanks a lot for your help ! your range idea it's pretty clever, it will already be much better than checking by hand!

Thanks @MarineBergot please remember to "accept" if it's useful. You can always switch "accepted" posts if there's a better way posted later.

@M__ by the way, do you have any idea why my current research is not returning anything?

@MarineBergot I've answered in the comments of comments

@M__ thanks for your answer. I arrived in the same point but i need both, latest and refseq. For me latest is working with organism and no partial but refseq is not working even alone. do you have results with Microbacterium[Organism] AND refseq[Filter]? could it be related to version ?

@MarineBergot its the database. The output of the above search using the nucleotide database is here

source: Stack Exchange

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

sort: votes


text:
calculate the backchain RMSD between two pdb files by pymol in Python
i have two protein pdb files and want to calculate the backchain RMSD between them. As far as I know the GUI of Pymol can use align to calculate the RMSD (but I don't know if it is the backchain RMSD).I want to use the python code to use pymol to calculate the backchain RMSD. I know how to align them but don't know how to calculate the RMSD.
import pymol2
with pymol2.PYMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar')

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

The backbone can be specified using a selector for the names of these, which are N, CA (Cα) and C (C-prime).
out = pymol.cmd.align('foo and name N+CA+C and polymer', 'bar and name N+CA+C and polymer')
rmsd, n_atoms, n_cyles, n_rmsd_pre, n_atom_pre, score, n_res = out

If you do not wish the proteins to be moved you can use rms_cur, which is just the rmsd.

Thank you for your answer! It works well! however,i want to further learn related knowledge about pymol like this. Can you provide some web/cook book name?

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


text:
Arlequin - How to assign samples to populations?
I have converted a VCF file to the arlequin file format. It consists of 20 samples that should belong to three populations, along with a few thousands SNPs. However, when I load the arlequin file into the software, it shows I only have 1 population. Any way to fix this?
Thanks.

Can you please show a few lines from your Arlequin file?

Arlequin format is historically extremely difficult to manually manage. The group leader eventually wrote software to allow easy Arlequin format management and conversion. This is achieve by Population Genetics Data (format) Spider, or PGDSpider by Lisher and Excoffier (2018). All information is available here and the manual is here. I would personally avoid attempting to manually rectify population designation.


text:
Nextflow tuple from CSV cols
I am trying to create a channel with a CSV file and use it as tuple input. Yet the output I obtain does not contain the entire file but only the first element of the file (1st row).
the csv
dig,val
1,One
2,Two
3,Three
4,Four
5,Five

output :
cat digOne.txt
1
instead of the whole collection
digOne, digTwo...

 size_amp_ch = Channel
        .fromPath( params.file)
        .splitCsv(header:true).map{
            row  -> tuple(row.dig,row.val)
        }

process{
    input:

    tuple val(digits), val(values)

    output:

    path "dig_${values}.txt"

    script:

    """
    echo $digits >> dig_${values}.txt

}

thank you

I think the problem is just that you have extra whitespace at the end of the first line in your CSV file. Either remove the whitespace manually, or use the strip parameter with the splitCsv operator to remove leading and trailing whitespace from values. For example:
params.file = './file.csv'

process test_proc {

    input:
    tuple val(digits), val(values)

    output:
    path "dig_${values}.txt"

    """
    echo "${digits}" > "dig_${values}.txt"
    """
}

workflow {

    Channel.fromPath( params.file ) \
        | splitCsv( header:true, strip:true ) \
        | map { row ->
            tuple( row.dig, row.val )
        } \
        | test_proc \
        | view()

}

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [stupefied_boyd] DSL2 - revision: 2a9d8dc414
executor >  local (5)
[4f/a73196] process > test_proc (2) [100%] 5 of 5 ✔
/path/to/work/9d/4905d0236788c1a3733da5358b7fc4/dig_Four.txt
/path/to/work/f9/d04e8a1c2265d73cde626703317df7/dig_Five.txt
/path/to/work/ff/522238f87a5641fdc7adce3e4d1cea/dig_Three.txt
/path/to/work/e9/9b7ec2a65d45986fd8944da4a2054c/dig_One.txt
/path/to/work/4f/a73196fbd56db3c9647fa7b072edfd/dig_Two.txt

Hey thank you so much for the answer, it seems that your answer is not solving my problem, I wonder if my test file is causing this. thank you so much for the help

@LoïcBorcard So what is the problem exactly? What are you trying to achieve or what are you expecting?


text:
Energy minimization of molecules
I have 3000 molecules which I need for docking but before that I need to perform energy minimization for all the molecules and it is taking a long time because I am doing it one by one. So is there a way that I can select multiple molecules and do the energy minimization

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Could you kindly consider improving this question? For example, what are the "molecules" in question, how is energy minimisation being performed and possibly an example. The answer is scripting via a loop and there are a number of members here with the skills who could help, but further information is likely to be needed.

As @M__ has mentioned, the terminology in your Q is a bit confusing.
A molecule can be the small molecule (ligand, compound etc.) or a polypeptide chain (protein).
Conformer generation of a compound is commonly not called 'energy minimisation' for technical reasons. For example, you normally generate in rdkit a conformer via AllChem.EmbedMolecule(mol), but you could do it via
p = AllChem.MMFFGetMoleculeProperties(mol, 'MMFF94')
ff = AllChem.MMFFGetMoleculeForceField(mol, p)
m = ff.Minimize()
assert m == 0, 'failure'

The latter would make sense if you want to add constraints (e.g. ff.MMFFAddPositionConstraint(0, 0.1, 10)).
There is also DFT (e.g. Gaussian or Schrödinger Jaguar), wherein the conformer is generated from a forcefield, but using a slimmed-down quantum theory to described orbitals as opposed to atoms. This is much slower than the above which use atoms as rigid particles tethered together (molecular mechanics).
The EmbedMolecule uses ETKDG method, which is a stochastic search method that cheats by using knowledge derived from crystal structures (cf. CCDC), so is way faster. It will be based on crystal structures, where the latter makes the compound best suit the terms of the MMFF94 forcefield.
If you are using a conformer generator, or similar, then the only quicker method is to retrieve the crystallographical structure.
If you are talking about minimising a whole protein, every forcefield will give slightly different results and there is no DB you could use. However, if it is for a snapshot/pose calculation, a common hack is do a local minimisation around your target and use implicit solvent which will be much faster. You have added the tag MD, which would suggest FEP calculations, for this local minimisation would not work.


text:
How much does Nanopore cDNA Sequencing Cost?
[this question is based on a question that was asked on Reddit]
We're interested in doing whole-transcriptome cDNA sequencing of 24 mouse cell lines, and are deciding between Illumina and Nanopore sequencing. What are the costs involved in Nanopore sequencing, and how much hassle is it to analyse Nanopore data?

The standard PCR-cDNA barcoding kit is \$750 USD, and has enough reagents for six runs, with each run using up to 24 samples (i.e. 144 samples total). Each run also needs one MinION flow cell (\$900 USD).
The total cost including external reagents is naturally dependent on the number of samples, and gets cheaper as the number of samples increases. For 24 samples, it's about \$5200 NZD (\$3100 USD). See below for more details.
If you could outsource this to a lab that has a PromethION (or P2 Solo), the sequencing would be a great lot cheaper. PromethION flow cells start at \$3600 for a pack of 4 (i.e. the same unit cost as MinION flow cells), but have about ten times the yield of MinION flow cells (which means more samples, or greater depth).
From a bioinformatics perspective, there's more work to do with nanopore reads because there are no commonly used / established pipelines (so I've had to create my own), but the results are much cleaner once they've been processed.
Cost Breakdown
For cDNA sequencing in Aotearoa, MinION works out cheaper than Illumina when compared against a service centre, about two thirds the price for us. Here is the reagent list I wrote up a couple of months ago for sequencing from 30 whole transcriptome samples, mouse or human (prices in NZD). This is with an ideal read threshold of 1M reads per sample (which I have found to give good quantitative results for differential expression). You could probably get away with lower yield for microbial sequencing:
ONT reagents:

1 x PCR-cDNA Barcoding Kit 24 \$750 USD [\$1210 NZD]
1 x Rapid Adapter Auxiliary V14 \$288 USD [\$481 NZD]
1 x R10.4.1 flow cell @ \$900 USD [\$1503 NZD]

Total ONT cost: \$3194
Other reagents:

10 mM dNTP NEB N0447 [\$145 NZD]
Maxima H Minus RT ThermoFisher EP0751 [\$346 NZD]
RNAclean or AmPure XP beads [AmPure XP: \$800 NZD]
Exonuclease 1 NEB M0293 [\$160 NZD]
RNAseOUT ThermoFisher 10777019 [\$569 NZD]

Total external reagent cost (excludes nuclease-free water, tubes, pipette tips, gloves, ethanol, DNA QC): \$2020 NZD
[As well as whatever you need to extract 200ng total RNA per sample]
Our current protocols for this PCR-cDNA barcoding + V14 RA approach can be found here.
The comparative Illumina cost according to the service centre we asked a couple of months ago was \$185 NZD per sample + \$5000 NZD per flow cell (or $9440 NZDfor 24 samples).
For almost all of these external reagents the amount ordered is an absolute minimum (from memory I think most are something like 50 reactions), so doing a second equivalent experiment would mostly only need the flow cells to be purchased.
In addition to the cost, turnaround time with MinION is quicker (for us it's a couple of weeks versus a couple of months for Illumina [with most of that time spent waiting in a queue]), and the results are nicer to analyse: full length transcript sequences with very even coverage profiles.
I'm aware that sequencing is more expensive in Aotearoa than it is in other countries, so the cost / time difference may not be as extreme elsewhere.
Bulk discounts for ONT flow cells
The prices above are for ordering reagents and flow cells one at a time, using the most expensive reagents. I call this a "grant funding" budget, because it's the quote I would give to people who are looking to fund research. There are ways to reduce costs when organising sequencing in bulk, and I've already hinted above at the extra capacity in reagent kits.
However, the biggest savings for nanopore sequencing can be achieved by ordering lots of flow cells at once. For example, ordering 48 MinION flow cells means that the per-flowcell price is \$500 USD instead of \$900 USD. That's 48 flow cells spread over the course of an entire year, with costs paid on delivery. If you can encourage multiple research groups to put in a combined bulk order, it's not too far out of the realm of possibility.
The long-read advantage
I've been told that a researcher is able to get 50M paired-end Illumina reads from a sequencing service centre for about $100 CAD.
50M paired end reads is about 10 Gb, which can be achieved on a MinION run (\$600 - \$1000), or about 1/10th of a PromethION run (\$850 - \$1500).
Even given those prices, it doesn't take much of a long-read advantage for MinION to become the most economical option. My own studies suggest that 1M full-length reads from nanopore sequencing (i.e. about 1/10th of a MinION R9.4.1 flow cell) gives results that are at least as good as 20M short paired-end Illumina reads.


text:
Is there a way to calculate the number of parsimony informative characters in a morphological dataset in R?
I have a dataset of binary morphological characters and would like to know how many of them are parsimony informative. Is there a way to calculate this with an R package?

The ape library (Analyses of Phylogenetics and Evolution) using the MPR method (Most Parsimonious Reconstruction).
Information on ape is here.
The format
install.packages(“ape”)
library(ape)
MPR(x, phy, outgroup)

x :   a vector of integers.

phy : an object of class "phylo"; the tree must be unrooted and fully dichotomous.

outgroup : an integer or a character string giving the tip of phy used as outgroup.

Rooting the tree is recommended. In the example rooting isn't used but an outgroup is used as taxa t1, that will be fine.
Generating random data,
vectx <- rpois(50, 1) # make random data
tree <- rtree(50, rooted = FALSE)  # also rtree(50) make random tree
MPR(vectx, tree, "t1")

Thats it, there's nothing more. Data, tree = steps.
Clarification the question is number of parsimonious characters. As a strict definition the number of informative characters is the same as the number of steps and thats what MPR will describe ... thats easy.
There are other possible meanings

The number of parsimonious characters/steps to make a tree is calculated via MPR as above.
if its the actual changes at a given position in the alignment (or vector) for each node in the tree - thats the answer below - but thats complicated.
If its simply is this site (not character) is parsimoniously informative and this one isn't then the "alignment" or vector could be used one character position at a time and subject to MPR and the step length would determine how parsimoniously informative each site was ... thus point 3 is easily doable in R

Note Calculating the actual sites on the alignment and associated changes per node is easily done outside R, but within R I don't know if thats doable.
Wait, I'm wrong it's doable but I would envisage its difficult to implement. Calculating parsimonious characters for each given node in a tree - yes I do that - doing that in R ... errhhh thats not me. There are Pythonic ways in, but I've not implemented them.


text:
How do I measure the proportion of different bacteria in a sample from a high-throughput sequencer?
[this question is based on a question that was asked on Reddit]
I have sequenced some [mostly cell-free] DNA from a sputum sample on a nanopore sequencing machine, and would like to know what is in it (i.e. the proportion of different bacteria). The run produced 2 million reads with an average read length of about 200bp.
What's the best way to go about doing a metagenomic classification of sequenced reads?

If you've got a fast solid state drive with >200G of space, then kraken2 + bracken works really well on both long and short reads (e.g. Nanopore, Illumina). It's also specifically designed for metagenomic analysis, unlike other general aligners like STAR, bwa, bowtie2, etc., so you get tree-based classification which can deal with partial / incomplete matches.
First, use kraken2 to map reads to taxa in memory-mapped mode (to reduce system memory consumption):
$ kraken2 --threads 10 --db /mnt/ultra_fast/kraken/metagenome --quick --memory-mapping reads/*.fastq.gz --report kraken2.report.txt > kraken2.out.txt
Then, use bracken to adjust the counts based on mapping / library probabilities:
$ ~/install/bracken/bracken -d /mnt/ultra_fast/kraken/metagenome -i kraken2.report.txt -o bracken.report.txt
I use pavian to look at the results, because it creates quite understandable Sankey plots:

But inspecting the _bracken_species.txt file is a pretty good text-based way to do it. If there's substantial amounts of a particular bacteria in the sample, it should come out in the results:
100.00  232939  0       R       1       root
99.99   232912  0       R1      131567    cellular organisms
68.52   159618  0       D       2           Bacteria
67.67   157635  0       P       1224          Proteobacteria
67.47   157152  0       C       1236            Gammaproteobacteria
61.82   144001  0       O       91347             Enterobacterales
61.79   143941  0       F       543                 Enterobacteriaceae
32.28   75183   0       G       561                   Escherichia
32.26   75147   75147   S       562                     Escherichia coli
0.02    35      35      S       1499973                 Escherichia marmotae
0.00    1       1       S       564                     Escherichia fergusonii

kraken2/bracken indexes can be downloaded from here (the build process takes ages; it makes more sense to use the pre-built indexes):
https://benlangmead.github.io/aws-indexes/k2
As far as I'm aware, kraken2 can't do strain-level assignment with the most commonly-used databases (although that may depend on what information that is stored in the database). For strain-level assignment, it seems like MetaMaps may work:
https://doi.org/10.1038/s41467-019-10934-2

Very nice, not possible to classify at strain resolution I guess?

Not as far as I'm aware; you could try metamaps, but I have no experience with that.

If you have bracken data from multiple samples and you’d like to compare them with relative abundance stacked bar plots, I made a shiny app for that: https://github.com/acvill/bracken_plot

Also, I believe taxonomic level S1 represents sub-species (strain) classification. Here are some example bracken output files representing different taxonomic levels, includingS1.

Cool question and answer ... thanks. Also thanks to @acvil


text:
Can I use palantir for trajectory analysis in a a merged object?
I was asked to use palantir for trajectory analysis. I have 4 samples that I have preprocessed using scanpy. Can I run the palantir analysis on the merged object or do I need to unmerge them?
Thank you

Thats a good question (but a bit brief). Palantir is here.

I have found the palantir code and I am trying to use it for my data. I was just asking if it makes sense to use it on the merged object of the 4 samples or if I should split them.


text:
How to make a UMAP for single cell data and color cells by average expression of a list of genes in scanpy?
I would like to make a UMAP where the cells are colored by the average expression of the bulk signature genes but I am not confident that I did it correctly. I would like to use scanpy for it.
I did the below:
bulk_de_genes_up_list = bulk_de_genes['Gene'].tolist()
#Subset the data based on the list of genes
adata2 = adata[:, adata.var_names.isin(bulk_de_genes_up_list)]
average_expression = adata2.X.mean(axis=1)
adata2.var['bulk_de_gene_average'] = average_expression
sc.pl.umap(adata2, color='bulk_de_gene_average', cmap='viridis')

I do get a UMAP as an output but I am not sure if it is done correctly. I am mainly worried about average_expression = adata2.X.mean(axis=1)
Is that the correct way of calculating the mean of the gene expression per cell?
Thank you


text:
Get a comprehensive list of virulence factors for certain or all pathogens?
Is there a database with all known proteins/genes of specific pathogens such as Staph aureus and annotated whether it is a known virulence or resistance factor that idealy can be downloaded in table format?

Please consider accepting MP's answer. For my own stuff I found it really useful, its a good a question and really good answer.

I suggest that you are looking for a resource like VFDB (virulence factor database). That table does annotate by bacterial species, so you can filter as you like.
You can download a description table or FASTA sequences from their downloads page (see screenshot). Last update was 4 days ago, so they should be quite current.

For resistance genes, you could use something like CARD (Comprehensive Antibiotic Resistance Database), which also has a downloads page where you can download a similar TSV table or FASTA sequences.
Both resources include web tools for predicting the presence of the elements of interest in sequences.
If you are looking for a state of the art prediction tool for such elements in genomes, you might look into seqscreen from Todd Treangen's lab.
For a general purpose microbiological annotation and analysis resource, you might also look into BV-BRG, which succeeds PATRIC. I have only used PATRIC but BV-BRG seems very thorough.
I think that you are unlikely to find a resource that has non-virulent/non-resistance genes annotated.

Wow what a cool answer. Got to be the accepted answer, but lets wait. There are a number of microbiological bioinformaticians here.

Thanks for the input. I saw the VFDB database, but I had concers how up to date it is. I am going to check out the other resources you mentioned.

@Soren added a note- VFDB was last updated 4 days ago. See screenshot.


text:
Primer trimming-fasta files
My goal:
I want to trim off the primers (Forward : CGAGAAGACCCTRTGRAGCT, Reverse : GTTGGGGYGACCNYGG) from a fasta file with a lot of dna sequences allowing for some (e.g. 3) mismatches (identity).
Input sequences
Seq1:

CGAGAAGACCCTATGGAGCTTAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTGGTTGGGGCGACCGCGGAACGT

Seq2:

ATGGCCCATCCCTCACAGCTAGGATTTCAAGATGCAGCTTCCCCAGTTATAGAAGAACTTCTCCACTTTCACGACCATGCCCTAATAATCGTTTTTTTAATTAGTACACTAGTA

Seq3:

CATAAGACGAAAAGACCCTATGGAGCTTTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTGGTTCGGGCGACCGCGG

Output sequences
Seq1:

TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG

Seq3:

TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

Seq2 is absent because it has no match for the primers so it should be discarded.
In seq 3 there is a mismatch in Forward (CGAAAA...) and in Reverse (GTTCGG...).  I don't need the primers to have 100% identity.  85-90% is OK.
Issues:

primer sequences can be anywhere in each sequence, not just in the beginning.

The sequence between the primers is not constant.  it is approximately 80 bases but not exactly 80.

I don't know how to put the primers' identity parameter.

How can I search for R or N nucleotides (IUPAC)? Because I don't know how, I wrote the primers with specific nucleotides and not R or N.

This is the code I have so far, but it has the problems mentioned above:
from Bio import SeqIO

def trim_adaptors(records, adaptorF,adaptorR):

    len_adaptorF = len(adaptorF)
    len_adaptorR = len(adaptorR)

    for record in records:
        indexF = record.seq.find(adaptorF)
        indexR = record.seq.find(adaptorR)
        if indexR == -1  or indexF == -1:
            # adaptor not found, so won't trim
            print(record.id, "no primer/s match")
            continue
        else:
            # trim off the adaptor
            yield record[indexF + len_adaptorF:indexR + 1]

original_reads = SeqIO.parse("trimming_testfas.fas", "fasta")

trimmed_reads = trim_adaptors(original_reads,"CGAGAAGACCCTATGGAGCT" ,"GTTGGGGCGACCGCGG")

count = SeqIO.write(trimmed_reads, "trimmed.fasta", "fasta")

print("Saved %i reads" % count)

This doesn't have to be a script in Python, I am also open to using existing tools.
Thank you a lot in advance, hope you got me :)

The IUPAC ambiguity codes can be thought of as regular expression character classes. R matches any purine, so [AG], Y matches any pyrimidine, so [CT], and N matches anything at all, so . or [ACTG]. So those are easy to handle: instead of  CGAGAAGACCCTRTGRAGCT you can use the regular expression CGAGAAGACCCT[AG]TG[AG]AGCT, and instead of GTTGGGGYGACCNYGG you can use GTTGGGG[CT]GACC.[CT]GG.
Allowing for random mismatches is trickier though and since you're open to using existing tools, I strongly urge you to use something like cutadapt instead of trying to reinvent the wheel. So, install cutadapt on your system (see here for instructions) and then run:
cutadapt -a "CGAGAAGACCCTRTGRAGCT...GTTGGGGYGACCNYGG;max_error_rate=0.15;" --discard-untrimmed -o output.fa file.fa

That assumes your input sequences are in fasta format in the file file.fa and  will create the output file output.fa:
$ cat output.fa
>seq1
TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG
>seq3
TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

The options used are -a $primer1..$primer2 which tells cutudapt these primers are linked adapters (https://cutadapt.readthedocs.io/en/stable/guide.html#linked-adapters) and then we also add max_error_rate=0.15; to allow up to 85% identity matches. Cutadapt already knows how to handle IUPAC ambiguity codes so nothing special is needed for those.
I am not sure how this will handle cases where only one of the two sequences is present, you'll need to test for that. I expect that they will be discarded since we're specifying linked primers but I would double check to be sure. In any case, cutadapt is a very powerful tool, so I'm sure there will be a way to handle that in the way you want if you check the documentation.

Thank you a lot! Sequences with only 1 primer (and the 2nd is absent or more mismatches than the given threshold) aren't discarded.  cutadapt cuts the one primer which is present and let the sequence in the dataset which i dont want to.  But generally it works.  Σε ευχαριστώ θερμά! @terdon


text:
Remove files after process is terminated
In my pipeline nextflow, I create a channel with channel.fromSRA and the channel contains lots of heavy files.fastq.gz. Then I have a first process to unzip files and transform them in files.fasta and a second process to count kmers. Is it possible to remove files.gz after first process and remove files.fasta after the second process. Can someone explain me how to do? They are too heavy to keep them in the work directory. Thanks.
In my pipeline I use kmer-ssr to count kmers present on my genome sample.
First I create a reads_ch = Channel.fromSRA
In this channel there are 12 files (. fastq.gz) of 20 Gb each one.
Then the first process unzips them and transforms them in files.fasta.
Now in my work directory there are 12 more files of 80 Gb each one. So the work directory is now 1200 Gb.
The second process counts kmer present in reads and creates 12 files of 100 Gb each one. The work directory will be 2200 Gb.
The work directory is too heavy. I want to know how I can suppress the current file. fastq.gz after first process

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Certainly this will be doable. Could you provide the code? Or else sufficient amounts of code to describe the nested loop? I will provide a link below as typical example of a NextFlow question.

So this is yesterdays question here as an example of a standard NextFlow question on the site.

Look into the "temporary" attribute of the inputs that should help you in the right direction.

Unfortunately, this remains an open issue with Nextflow: https://github.com/nextflow-io/nextflow/issues/452
However, writing large amounts of uncompressed data to disk is usually best avoided. If it can't be (e.g. maybe the file needs to be read multiple times), you'd need to merge the two processes so that the writing of the temporary file and the running of the tool is done in a single process. This would then let you remove the temporary input file as soon as it's no longer required.
Another way would be to pipe the uncompressed and FASTA formatted reads directly into your tool. Not all tools support this functionality obviously, but Kmer-SSR can indeed receive input from STDIN, which completely obviates the need to write any temporary files. Note that by default, Kmer-SSR writes also to STDOUT and this lets you write a compressed output file:
process kmer_ssr {

    cpus 4

    input:
    tuple val(sample), path(fastq)

    output:
    tuple val(sample), path("${sample}.tsv.gz")

    """
    zcat -f "${fastq}" \\
        | sed -n '1~4s/^@/>/p; 2~4p' \\
        | kmer-ssr -t ${task.cpus} \\
        | gzip \\
        > "${sample}.tsv.gz"
    """
}

Since the above uses a pipeline, make sure to add the pipefail to your command. One way is to simply apply it to all processes by adding the following to your nextflow.config:
process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}

From the bash man page:

pipefail
If set, the return value of a pipeline is the value of the last
(rightmost) command to exit with a non-zero status, or zero if all
commands in the pipeline exit successfully. This option is disabled by
default.


text:
How can I detect MCPyV in a FASTQ file?
I have received two NGS files from an NGS company, both are FASTQ files that correspond to reads from a tumor sample. I heavily suspect that MCPyV is present, and I am hoping to identify it.
What downloadable software can I use for viral-host integration detection? I have access to both Windows and Mac OS, but have no experience with Linux.

You could start by trying out a metagenomics analysis program like kraken2. There are a few databases that combine bacteria, fungi, and viral searches.
One benefit of using kraken2 is that it has options to output the most likely taxa of matches at the sub-read level, which can help with identifying integrated or chimeric sequences.
I've written up the higher-level approach in more detail in this answer:
How do I measure the proportion of different bacteria in a sample from a high-throughput sequencer?
The sequence-level hypotheses about kmer origin can be found in the direct output of kraken2 (last column):
$ kraken2 --threads 10 --db /media/minion/minion_temp/db/kraken/metagenome --memory-mapping fastq_pass/barcode0?/PAK83314_pass_barcode0?_*.fastq.gz --report kraken2.report.txt > kraken2.out.txt
$ head kraken2.out.txt
C       a3e0b7c9-e6f3-40bf-900c-b1c59eb4977d    9606    208     0:70 623:3 0:53 9606:1 0:6 9606:32 131567:5 2759:4
C       8030d33f-3838-4ed3-bab6-5f3b79fd059e    575614  164     0:68 623:1 0:46 575614:2 0:13
C       5d87a648-799d-451f-b747-32a0fdf2517a    9606    190     0:108 9606:5 0:43
C       8aba667d-4a5b-4f8b-9d18-1b5f761b87b2    9606    267     0:124 9606:9 0:99 1:1
C       1b3501c2-6766-4ee8-b984-9e92470f37da    9606    335     0:70 28901:5 0:12 35885:1 0:23 9606:16 0:83 9606:17 0:6 9606:2 0:9 9606:5 0:11 9606:41
U       75c2625d-ab87-4614-8d0d-c9d12c7d1e19    0       242     0:208
U       902b0ec1-80d9-4e07-80c8-29a79191d931    0       171     0:137
C       56b7c8a8-4cfd-44ec-8922-62871b29d31c    9606    338     0:103 9606:5 0:22 9606:5 0:1 9606:13 0:55 9606:27 0:38 9606:3 0:32
C       68bf6718-2692-46e7-ba2d-f45952fd6138    9606    197     0:68 1:9 0:38 9606:19 0:29
C       81f0f4ae-3bb5-4b94-89ff-152ff3b6f5d7    9606    238     0:146 9606:46 0:12

I'm aware that you're interested in looking for integration of a single viral genome. This can be specifically looked at with kraken2 because it can handle custom-built indexes as well. However, the benefit of using a comprehensive index (e.g. the MinusB database from here) is that it will make it much easier to exclude the other numerous endogenous viruses that are already integrated into the genome, whereas using a MCPyV-only index may lead to false-positive matches to sequences that are from another virus.

FYI (non-OPs) it is very likely this sample has already been diagnosed and proven positive via a PCR test. I don't know the OP, but I am pretty familiar with oncovirus studies. To be clear, my guess is that PCR positive - doesn't mean cancer.

What I think is occurring (if I'm wrong hopefully the OP will say) - is the  OP has a cancer patient(s) and has good evidence that the underlying cause is MCPyV: which is a very specific oncovirus.
What they want to do is identify viral integration in the patient genome as evidence of oncogenesis. To put this general field in perspective I think this issue is literally the rise and fall of many investigators and companies in cancer (again not focusing on a singular virus per se but generically on the oncoviruses). Keep in mind I have an explanatory bias towards viruses.
The molecular genetics of integration is likely to be  complex.

In summary, the OP will particularly want to do (I suspect) is map the viral integration site. What is difficult is not diagnosis of the virus because at any time post-any sort of approximate assembly the genome can be locally blasted against the reference virus and its likely to have already been PCR'd. It's getting an assembly good enough to understanding the integration site - to demonstrate integration has occurred and thats is likely difficult (but I hope I'm wrong). Inverse PCR was used for this sort of thing - but thats a wetlab approach and PCR-based strategies with traditional Sanger sequencing can still be very useful here, but its laborious but was/is a  traditional approach.
Thus, I suspect the OP is from a classical virological background venturing into NGS. I noted the OP used a human reference genome hg38 in a previous post.
If thats correct then ... the big issue in my opinion is obtaining a correct assembly of the putative virus-human genome O Steps would be,

QC
Assembly e.g. via bwa-mem2 ... O O
Construct a local makeblastdb database of the genome then ..
blastn between the human genome and a viral query to identify the integration site
Perform something like a gatk analysis to have a different perspective on the integration site via VCF
hone in on the contig via IGV for optical analysis and study how the virus has integrated

O My eukaryotic genomics is outdated and the above approach is summarised in Genomics in the Cloud by Geraldine A. Van der Auwera, Brian D. O'Connor published through O'Reilly Media. I just focus on small genomes.
O O and NOTE I cannot remember if bwa-mem2 uses a reference genome e.g. hg38 type thing (Auwera and O'Connor will describe it). If it does the assembly of the integration site, in my opinion, would be compromised and alternative assembly strategies might (probably would) be required. I don't think this will be trivial - I think it would need a skilled eukaryotic bioinformatician who could be told about the biology of viral integration. Keep in mind, eukaryotic genome assembly is not my thing.
However, riches await the investigator if this can be satisfactorily solved IMO.

To answer @gringers point, Kraken2 would be useful for viral integration per se  and yeah the results would be interesting because there's quite of them ... but in this case - I think - its an exact genetic sequence of a single species, thus I think direct blast-type blasting or minimap2 mapping is a better strategy than Kraken2. Just my opinion.


text:
DNASTAR viral-host integration assembly keeps failing
I have two NGS files from an NGS company corresponding to the sequencing data from a tumor sample as follows:

TB_7710391_R1.FASTQ.gz
TB_7710391_R2.FASTQ.gz

I have downloaded the genome for MCPyV as a file called “MCPYV_Sequence.fasta”
I have downloaded the Hg38 as a file called “hg_38.fa”
I have input the host genome as hg_38.fa, the viral genome as “ MCPYV_Sequence.fasta” and the sequence data as both FASTQ files.
I run the assembly in the cloud, and I receive back an error stating that the assembly as failed.
What am I doing wrong?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

This question needs more detail. What cloud service are you using? How are you using it? Why are you trying to do an assembly? What results are you looking for?


text:
How do I filter reads from a bamfile based on 5' and 3' ends?
I have data in paired-end as well as single-end format for the same sample.I would like to split those two files into two per each, where one contains sequences that have 5' end and the other one has 3' ends.
So the paired end data has R1 and R2, if R1 is first, it has 5' ends and if not, it has 3' ends? And I'd do the same thing for the single end data.
I'm using samtools and pysam.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

It's unclear what you're asking. Sequences are almost always written out to files in the 5' to 3' direction.


text:
Download Addgene database
Is anyone aware of a way to download all vector sequences from addgene at once? As a non-profit, I'd assume they shouldn't be opposed to making this data available, but it doesn't seem to be possible. Might anyone have scraped it at some point?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

That is a very good question for Addgene. You could contact them directly.


text:
How to calculate cell type percentage in every sample
I have a Seurat object (metadata) with the single R samples consisting of cell types and sample types columns. I am trying to make a table that has a sample and percentage of cell types for each sample so that I can run a t-test across the samples that came out of 3 prime vs multi-omics.
I have a seuratobject meta data with singler_labels (cell labels like chondrocytes, astrocytes, macrophage...), sample names (S0030, S0031, S0032.....), cell id/barcodes, and everything a seurat object has. I want to make a table with samples and percent of cell types per sample just to see what sample has what percent types of cells.
How do I do this?
Thank you!

We can't really help you out without looking at your Seurat object. Can you show us what all the metadata in the object looks like?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Hi, do you have replicates for each condition you want to test? Based on your short description I assume you have but it would be useful to specifically mention this. Be careful when comparing such data because proportions are extremely variable from run to run

If your cell types and sample names are in separate metadata variables attached to the Seurat object, then you can use table to count up the pairings:
table(pbmc.seurat$singler_labels, pbmc.seurat$sample_name)

This will output a matrix containing these counts

The following package performs this type of analysis and can be directly used on a Seurat object:
paper: propeller: testing for differences in cell type proportions in single cell data
github: Speckle
You mention that you would like to do some statistics so keep in mind that this requires replicated measurements to handle the large variability in proportions between samples. Your metadata should contain:
celltype, sample, group_to_test
once the package is installed you can simply run the propeller function and set group to technology (in your post you mention you want to evaluate 3' vs multi-omics):
propel <- propeller(clusters = object@meta.data$celltype,
                    sample = object@meta.data$sample,
                    group = object@meta.data$technology)

This will return the statistics: t-test for two groups or ANOVA for multiple testing.
The following function will return the transformed proportions that can then be used in limma for example to test for more complex designs or simple sample-wise plotting using ggplot. For transform you will offcourse have to decide which transformation is applicable to your data.
getTransformedProps(clusters = object@meta.data$celltype,
                             sample = object@meta.data$sample, transform="asin")

There are a variety of other packages that allow for more complex analysis and a recent benchmark can be found here


text:
Why might parabricks not be writing output to the output file I provided?
This feels like a silly question, but I am testing parabricks bammetrics and nothing is being written to my output file. This is the command I am using:
docker run \
  --rm \
  --volume ~/workdir:/workdir \
  --volume ~/ouputdir:/outputdir \
  -w /workdir nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 \
  pbrun \
  bammetrics \
  --ref /workdir/reference.fasta \
  --bam /workdir/cram.cram \
  --interval-file /workdir/regions.bed \
  --out-metrics-file /outputdir/metrics.txt

But metrics.txt doesn't have anything written in it when I open it after the command is run. Am I doing something wrong? I am very new to parabricks, so I assume I am missing something simple
I apologize if this is the wrong place to ask this question.

EDIT Summary
For future users who encounter the problem the solution is,
docker run --platform linux/amd64 --rm -it -u 0 --entrypoint /bin/sh  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

I am unhappy with the practice of limiting Docker access via blocking standard runtime: there are statements of intent on SE as a strategy to attempt to restrict Docker access. Parabricks isn't an isolated incident of deploying this setup within bioinformatics. However, I checked their forum and there is no statement of intent that it's deliberate. Parabrick's copyrights are here. There would be issues if there was ever an official statement to this extent. They would likely argue it was for other reasons to automate stuff.

Original answer ...
I don't know parabricks, but I know Docker. What you are attempting is to run Docker from the outside ... I would get on the inside to view the error log. Whats happening is the error is not being passed outside the container. Getting on the inside requires
docker run  --rm -it --name mydocker nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

This should give you access to the container. If it doesn't post back and I'll forward a different command that will. It's impossible to block access to a Docker container albeit some people try.
For simplicity is worth trying
docker run  --rm -it --name mydocker clara-parabricks:4.0.1-1

Anyway, if you are using this approach you'll need,
docker cp *outside container* *insider container* nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 and vice versa

... this is an amazing command and allows you to move data in and out of the container. You've got to get the data on the inside of the docker and docker cp does that very well. I would then look to running the whole thing from the inside of the container. Thats exactly how I run Dockers, then I exit and return when it's done. There might be a bit of clean via this approach, I'm not sure --rm might take care of it.
You'll also be able to look at the programs and copy it out onto hard drive. This last bit may be not advisable - the whole point of Docker is to protect your system from their algorithms ... but you can definitely do it and if you get it paired up with NextFlow thats cool.
I don't particularly want to pull the Docker container because they're often huge and I've no interest in the algorithm albeit I've spotted its at nvcr.io/nvidia/clara/
Docker, to be honest it's lazy by the developers. Docker is - techie packaging. So it's good for reviewers and behind web firewalls. It's useful but like you've just experienced  if it flops you're stuck. The one positive thing about Docker (apart from DevOps) is it's generally tolerated to package someone else's algorithm as part of a pipeline and let a third party 'pull it'. This is because they'll be a Dockerfile which can reconstruct the whole thing without copyright infringement. So if that first developer objects ... the Docker developer would just issue the Dockerfile. conda and stuff max out on that concept.

There's a weird error as the OP described below.
Right Docker uses by default a tag latest that does not work for parabricks. The tag must be 4.0.1-1. Could you retry the above code please?

Update For good measure try the following:
docker run  --rm -it --name mystuff  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 --entrypoint /bin/bash

Might not be bash, thus
docker run  --rm -it --name mystuff  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 --entrypoint /bin/sh

They shouldn't throw an error, but you might still find you are outside the container after running it. Lets see.

Update
The error looks okay, because a runtime error is exactly what would be expected if they are attempting to prevent access. Therefore please try:
docker run --platform linux/amd64 --rm -it -u 0 --entrypoint /bin/sh  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

This, or a variant of it, should work.

Thank you for the response! When I ran that command it didn't work, I get this error:
edited
I've read that NVIDIA does not want to give access to the parabricks container

THANK YOU! I was able to get inside the container and run the command. For some reason, running it from within the container caused the output to be written to the metrics file... I'm not sure why it would be different, but I'm glad that it worked nonetheless. Thank you so much again for your help!

Just want to clarify that I read that nvidia doesn't want people to access the container somewhere on the dev forums, not in the documentation (and I was paraphrasing, the original dev response was "the best way to run parabricks is outside the container").  edited

Okay thanks @womy. Checked the forum, its fine.  Note some of the comments in this thread have been deleted or edited to comply with extended discussion SE stuff


text:
How to plot average gene expression in scanpy?
I would like to make a UMAP where the cells are colored by the average expression of the bulk signature genes but I am not confident that I did it correctly. I would like to use scanpy for it.
I did the below:
bulk_de_genes_up_list = bulk_de_genes['Gene'].tolist()
# Subset the data based on the list of genes
adata2 = adata[:, adata.var_names.isin(bulk_de_genes_up_list)]
average_expression = adata2.X.mean(axis=1)
adata2.obs['bulk_de_gene_average'] = average_expression
sc.pl.umap(adata2, color='bulk_de_gene_average', cmap='viridis')

I do get a UMAP as an output but I am not sure if it is done correctly. I am mainly worried about average_expression = adata2.X.mean(axis=1)
Is that the correct way of calculating the mean of the gene expression per cell?
Thank you


text:
Python list comprehension to calculate the set of categories for each row
My data = data
Model   A1BG    A1CF    A2M A2ML1   A3GALT2 A4GALT  foldchange
IM00499 N   L   L   G   N   N   0.285744
IM00499.1   L   N   D   G   D   D   0.005046
IM00980 N   N   N   N   N   N   -0.068949
IM01529 N   N   L   N   N   N   0.065914
IM01624 A   N   L   N   N   N   0.302100
IM01240 N   D   N   D   N   G   -0.428324
IM00807 N   L   A   N   D   N   0.257222
IM00999 N   N   N   N   N   D   0.056268
IM00983 N   N   N   N   N   N   0.107456
IM01076 N   N   D   A   N   N   0.218301
IM00803 N   D   N   N   N   D   -0.479368
IM00125 N   N   A   N   N   A   0.270428
IM00998 D   N   N   N   G   N   0.142748
IM00997 N   L   G   G   N   G   -0.154540
IM00996 N   L   N   N   N   L   -1.302023

The above data frame is divided into 5 categories = A, N, L, G and D
feature = ['A1BG', 'A1CF', 'A2M', 'A2ML1', 'A3GALT2', 'A4GALT']
Feature_level_partitions are arrays of combination of classification values
feature_level_partitions = ['A', 'N' , 'L',  'G',  'D'], ['D'], ['L', 'G']
I was using the below list comprehension to calculate Hedge's G.
        hedges_g = np.array(
            [
                Hedges_g(data[feature].isin(feature_level_partition), y)
                for feature_level_partition in feature_level_partitions
            ],
            dtype=np.float64,
        )

Hedges_g is defined as follows:
def Hedges_g(x, y):
    difference_of_means = abs(np.mean(y[x])[0] - np.mean(y[~x])[0])
    s_pooled_weighted = math.sqrt((((n1 - 1) * s1 ** 2) + ((n2 - 1) * s2 ** 2)) / (n1 + n2 - 2))
    return difference_of_means / s_pooled_weighted

I am interested in knowing how many models there are within each category. So for example if I divide up the data into ["D","L"] on one side, and the rest on the other, I want to know how many models have either "D" OR "L" and how many models have any of the other categories and then don't calculate Hegde's G if the number on either side is less than 2.
To answer the above query I tried the following code
[data[feature].isin(feature_level_partition).sum()
    for feature_level_partition in feature_level_partitions]

The above code gave me the number of models with the feature_level_partitions
However I also want to know the number of models without the feature_level_partitions
I also want to create a python function which calculates both above numbers in one line and fit this code as a if statement in my above list comprehension function

What you might think about @Megha is defining a new question in the context of both the comprehension and the models by defining a criteria for including 'foldchange' parameter.

What I'm not clear about is the definition of y - thats crucial to understanding what the code does. Just to note, although you've called the dataframe data, I've called it df (traditional shorthand). In terms of speed there will not be much in it between vectorised pandas and numpy. Explicitly using iterrows or else apply is when it slows down a lot.
Possibilities of the definition of y.
1.
y = df.drop(['Model', 'foldchange'], axis=1)

I think the above is more likely, or ...
2.
y = df['fold change']

What I think is occurring is np.mean is performing a boolean operation on the feature columns. What I didn't understand though was this ...
abs(np.mean(y[x])[0] - np.mean(y[~x])[0])

.. this is because for example model 'IM00980' albeit the mean is 0 this will always be 1 under the above equation. This is because y[~x] will be 1 for all values and thus abs(0 - 1) is 1. Equally if all columns in a row were 1 the output would also be 1 because abs(1-0) is 1. Anyway, what I would do is explicitly perform the operation within pandas rather than coercing np into a boolean operation. This is
import pandas as pd

df = pd.read_csv(r'/pathtofile/megha.csv', delimiter="\t")

feature = ['A1BG', 'A1CF', 'A2M', 'A2ML1', 'A3GALT2', 'A4GALT']
df=df.drop(['foldchange'],axis=1)
feature_level_partitions = ['L', 'G']
filtered = df[feature].isin(feature_level_partitions).eq(True).mul(1)
df2 = df['Model'].to_frame().join(filtered)
df2['mean'] = df2[feature].mean(axis=1)
print(df2)

Here ~x would simply 1 - mean as a new pandas column, i.e. df2['~mean'] = df2['mean'] - 1
The output of that code is:

-
Model
A1BG
A1CF
A2M
A2ML1
A3GALT2
A4GALT
mean.

0
IM00499
0
1.
1.
1.
0.
0.
0.500000

1
IM00499.1
1
0
0
1
0
0
0.333333

2
IM00980
0
0
0
0
0
0
0.000000

3
IM01529
0
0
1
0
0
0
0.166667

4
IM01624
0
0
1
0
0
0
0.166667

5
IM01240
0
0
0
0
0
1
0.166667

6
IM00807
0
1
0
0
0
0
0.166667

7
IM00999
0
0
0
0
0
0
0.000000

8
IM00983
0
0
0
0
0
0
0.000000

9
IM01076
0
0
0
0
0
0
0.000000

10
IM00803
0
0
0
0
0
0
0.000000

11
IM00125
0
0
0
0
0
0
0.000000

12
IM00998
0
0
0
0
1
0
0.166667

13
IM00997
0
1
1
1
0
1
0.666667

14
IM00996
0
1
0
0
0
1
0.333333

If this is correct it provides the basis to move forward for all other operations. Basically, what I'm saying is I think is good to break open the numpy into explicit pandas operations, what is happening is there's a mix of pandas and numpy and I would recommend going down one route or the other.


text:
Can we test the accuracy of Phred scores shown in FASTQ files
Recently I have ran a human WGS on the BGI DNBSEQ system, and their FASTQ quality scores seem to be quite impressive, where the Phred scores barely deteriorate along the read length when checked on FastQC.

Since the DNBSEQ technology is relatively new to me, I wonder if there's any bioinformatic approach to test whether the Phred scores listed in these FASTQ files are authentic or not? I'm kind of worried that their sequencer would just falsely assign high Phred scores to all bases.

How do the other FastQC metrics look? How does the alignment/mapping QC look?

Good question!  To Chris_Rand's point, you could try to test it by seeing how often a base call actually looks correct versus what each Q score implies, if you can gather that information separately.  My first thought was, what does BGI claim goes into their Q scores in the first place?  If you look it up for Illumina you can find at least general overviews of how they derive them from the image processing metrics like signal intensity and signal-to-noise ratio.  I can't find anything like that for DNBSEQ, so it seems like more of an open question there.

Yes, reported accuracy can be tested and compared with actual accuracy by mapping reads from a reference sample to a known-perfect reference sequence. It's usually done in aggregate, at a read level rather than at single-base level.
I like using Heng Li's gap-compressed identity statistic for this:

Alternatively, fish out a few reads at random, and use BLAST, which reports an identity score.
If the sample is not a reference sample with a precisely known sequence, then this approach will not work, because introduced variation will confound the comparison.
While it might be possible to create a new reference sequence from the reads, this introduces the possibility of systematic error in the sequencer that would be hidden by the new reference sequence, and when you're looking at error rates in the realm of q30 and above (as in your plot), those small differences are going to matter a lot.

Thank you very much for your answer! But how can the gap-compressed identity/BLAST approach distinguish between SNPs and sequencing error? If the gap-comprs identity is really low, should I compare the number of SNPs called with "normal" WGS sample? If in this case, the number of SNPs is also higher than "normal", then I can infer the low gap-compressed identity score originates from high number of SNPs in the sample? Or should I take a look at the allelic count in the VCF, since it should be small at the error base? Sorry for having so many additional questions, really appreciated your help!

Seems like you're asking a different question with that (i.e. "How can I trust variants discovered through high-throughput sequencing?"). The problems you state are definitely an issue, and they can't be resolved fully when using a non-reference sample as input.


text:
Should I compensate for two clusters of samples in microarray by hierarchical clustering and MDS?
I seem to have two groups of samples based on hierarchical clustering as well as MDS plot analysis using the sample label from the tree cut of the hierarchial clustering. I am analyzing microarray data. Should I adjust for the two different groups of samples in the model? For context, my covariate of interest is a continuous measure (air pollutant measure). This is a human cohort and not a case-control model and hence I was wondering if I should adjust for the different groups of samples.
[
MDS plot is here for the same samples:

Very interesting, if I understand correctly you are asking whether a straight hierarchical clustering should account for a continuous variable because it is resulting in discrete groups (clusters) which seems contradictory for a continuous variable?
If that is the question the answer is strictly no. Do not use weighting for this particular analysis (other regression style analysis yes - but not this one). Thus 'weights' are for regression analysis not for multivariate statistics. This clustering analysis is about assessing the correlation and not trying to impose it - thats regression.
The 'weighting' that can be used is different for a tree based analysis - it accounts for rate variation within the expression data - you are assuming rate homogeneity. That doesn't always work - but it's complicated and your data looks good as is.
Moreover, I would suggest your have 4 clusters. What you want to do is map the continuous trait across all branches 'tips' of your clustering. This will break up your continuous variable (hopefully) and looking at the MDS very much into discrete groups. That is the result, in fact it may not be 2 clusters but 4 clusters - map all the 'tips' against the continuous variable. There're statistically randomisation tree methods to assess the robustness of the observation - you can also use bootstrapping.
The question is whether these discrete groups are in fact the best description of the data. The MDS suggests that the continuous variable "air pollution" is in continuous within gene expression - using a combination of clustering and MDS you highlight this facet of the data. There alternative analyses namely regression or GLM to assess the relationship between the continuous variable and your results - the MDS is suggestive of that. However, it's not a 'and', 'or', 'right', 'wrong' relationship a continuous variable can be broken up into discrete groups: if there is a 1:1 relationship between the gene expression and the continuous variable "air pollution" (MDS analysis) - that is a better description, if that relationship doesn't hold breaking into discrete clusters is the best description thus far.
Essentially this is statistical modelling and the question is what is the best model to fit all the data. I'm personally not a fan of hierarchical clustering as a tree description of the data - but in your data it looks good.
At the heart of the question you are asking whether there are boundary effects in the data as described by the hierarchical clustering? Thus is there a certain level of pollution (or whatever) causing one gene expression boundary and then another and so on ... Perfectly possible - but central to that hypothesis is whether clustering is the best description of the data rather than a regression analysis, i.e. is it continuous or is it discrete.
Discrete vs. continuous  If it is discrete then k-means clustering would be a really, really good follow-up analysis to confirm the discrete "boundary" effects within the data. GLM is an analysis that might trap the data if it's continuous, there a lot of other methods in regression analysis. Basically if you map each "pollution" exposure against the MDS and if it's a perfect fit - trapping the result in a regression model outside eigenvalues would be the traditional thing to do. However, these days "unsupervised learning" (of which MDS is one method) is trendy so it may not be needed.
k-means clustering I would certainly consider k-means clustering in any case - this is because hierarchical clustering makes a rate assumption which might smash up the discrete groups. K-means can capture that without getting into complicated tree analysis. Basically, k-means you'll set 2 groups - assess the fit of the continuous variable, try 3 and then 4, how does the air-pollution variable fit? The data will also work well with supervised machine learning. However, it is genuinely continuous then regression analysis is best approach. With k-means you can also do parameterisation, i.e. whats the optimal number of groups represented by the data: if it comes back with '4' - which neatly segregate "air pollution" thats really good.
It's an exciting result. The way you write the result in the context of the conclusion is really important - i.e. the continuous verse discrete best fit analysis - in case you encounter a statistical modeller in peer-review.
My answer is a bit waffly, but is a representation of statistical modelling in this particular data, because two approaches may both be appropriate, but one is a better description than the other.

Note, if I've misunderstood the question and some sort of standardisation is needed between experiments - thats a different answer. Thus if you are saying that the groups represent different experiments, e.g. performed on different days or different arrays and there's risk of experimental artefact: thats about standardisation and the answer would be absolutely. Standardisation is essential to rule out experimental bias - particularly in microarray data.

Thank you for your very detailed answer. Highly appreciate it. I will slowly go through it and digest it. :)

Thanks let me know @siddharthadas


text:
Default placeholder for missing fields doesn't work with block in xtract tool
I am using the xtract tool of the Entrez Direct (EDirect) package of NBCI. I have a list of Accession Numbers:
U47804, U47803, U47802, U47801, U47800, U47799, U47798, X92938

and try to get information about:

Name of organism (Org-ref_taxname tag)
Product of gene (Prot-ref_name_E tag)
Public year (Date-std_year tag)

I executed the following command:
for i in {U47804 U47803 U47802 U47801 U47800 U47799 U47798 X92938}
do
  efetch -db "nucleotide" -id "$i" -format full; done |
    xtract -pattern Seq-entry -element Org-ref_taxname \
      -block Prot-ref_name -def "-" -element Prot-ref_name_E \
      -block Date-std -position first -element Date-std_year

The result I got:
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  1995

Because the X92938 doesn't have <Prot-ref_name>, so I was expecting to replace the missing fields with "-" by the -def tag of the extract tool. It means my desired result is as follows:
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1    -     1995

Can someone have experience with the NBCI tool may help me? Thanks a lot.

I saw this on stackoverflow but I couldn't figure out the problem. One of the great things about https://bioinformatics.stackexchange.com/ is that you can include heaps more detail here; if you could please provide more context it might help e.g. Why do you want this data? Have you looked at other sources e.g. the Rentrez package? Could you post-format the data using e.g. awk? Etc

Thanks for your helpful suggestion. According to your suggestion, I have used Rentrez and realized it is much easier to use and more powerful than the Entrez bash. So in some special cases like creating bioinformatics pipelines, using bash Entrez is a compulsory thing.

After trying myself in many ways, I have found that the method of creating a checkmark column when a desired column is recorded is more efficient and easier to understand.
Here is my bash script to solve my problem in a verbose way (I guess has others way shorter).
#!/usr/bin/bash

change_index(){
    local in="$(< /dev/stdin)" # Variable assignment from the stdin
    local col_num=$(head -n 1 <(echo "$in") |awk -F"\t" '{print -1+NF/2}')  # Number of fields

    for k in $(seq 0 $(($col_num-1))); do # Loop through number of field
        in=$(awk -v ia="$k" -F"\t"  '{sub("/[^"ia"]","/"ia"\t-\t/"((ia+1)),$(2*ia+1)) ;print $0}' <(echo "$in")) # Supplement missing fields
    done

    in=$(awk -v ia="$col_num" -F"\t"  '{sub("/[^"ia"]","/"ia"\t-",$(2*ia+1)) ;print $0}' <(echo "$in")) # Supplement missing the last fields
    in=$(sed 's/\/[[:digit:]]\+//g' <(echo "$in")) # Remove markers
    echo "$in"
}

for i in {U47804 U47803 U47802 U47801 U47800 U47799 U47798 X92938}; do \
    efetch -db "nucleotide" -id "$i" -format full; done |
    xtract -pattern Seq-entry \
        -lbl "/0" -element Org-ref_taxname \
        -block Prot-ref_name -lbl "/1" -element Prot-ref_name_E \
        -block Date-std -position first -lbl "/2" -element Date-std_year |
    change_index| column -t

Explanation
The main idea is that if a field is recorded there will be a marker field lay before the recorded field. This action is done by -lbl tag. And we have the following draft table:
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /2  1995

After having a draft table, and cleaning the table by change_index() function I wrote.


text:
How can I identify the virulence factors in my FASTQ files generated from nanopore sequencing bacteria DNA?
I have performed genomic sequencing using nanopore technology and upon analysing using BLAST in the command line I have identified the presence of E.coli. I was wondering if I can add some code to identify the known virulence genes of E.coli which are on the NCBI BLAST database? Thank you.

If you could kindly try to improve the question unless @pippo1980's answer is helpful. This is not a clear question because you appear to be asking a very specific question for undefined code for a database we don't seem to be aware of, i.e. a NCBI virulence factor database. It's not that it's not solvable - but without further information its really hard, e.g. present an example of the Blast result you want to explore.

you could try to submit your sequences to VFDB: a reference database for bacterial virulence factors ; paper here : reference paper.
use the VFanalyzer tool

VFanalyzer provides an automatic analysis pipeline for a systematic screen of known/potential VFs in given complete/draft bacterial genomes. Users can perform the analysis by one of the publicly available bacterial genomes either from our pre-downloaded genome list (Option 1) or via any valid GenBank accession(s) (Option 2). Or alternatively, users can also analyse any private genome data (raw FASTA sequence, predicted proteins, or annotated GenBank format genome) using Option 3. For query with raw FASTA genome, predicted CDSs will be identified using GLIMMER3 prior to the run of VFanalyzer.
Please note that, VFanalyzer only accepts query genomes from the genera of bacterial pathogens with full information in VFDB, and it's designed for complete or nearly complete draft genomes. So you are at your own risk if the query genome are derived from other bacteria or only partially sequenced.

NOTE: Due to connetion issue with NCBI, VFanalyzer can not fetch gbk file using Edirect correctly. Please download the file manually and use option3 to submit your jobs, we will find the other way to fix the problem ASAP.

Think you need to convert your FASTQ to FASTA first.
Other resource is : VirulenceFinder 2.0 or PathogenFinder 1.1;  paper ref : In Silico Genotyping of Escherichia coli Isolates for Extraintestinal Virulence Genes by Use of Whole-Genome Sequencing Data, PathogenFinder - Distinguishing Friend from Foe Using Bacterial Whole Genome Sequence Data


text:
How do I build the "Standard Kraken 2 Database"?
I am at this point in the GitHub tutorial: https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown#standard-kraken-2-database
It says to create the standard Kraken database, I use the following command:
kraken2-build --standard --db $DBNAME

It then says to change $DBNAME with the preferred database name/location. Is this something that I need to download on my own? Where can I find it, and how do I do it?
To provide context, I want to use Kraken2 to identify whether or not there is the presence of MCPyV in a tumor sample.

If you just want the viral database, read the 'custom databases' section (https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown#custom-databases)

\$DBNAME is a folder location on your computer; you specify it yourself, for example:
kraken2-build --standard --db /media/ip80/externalDrive/kraken_indexes/betterThanGenbank


text:
Kraken2 Standard Database failing to build (unexpected FTP path)
I am attempting to build just the standard Kraken2 Database by using the following command:
kraken2-build --standard --threads 24 --db $DBNAME

I am returned the following error a few minutes building it:
rsync_from_ncbi.pl: unexpected FTP path (new server?) for
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/762/265/GCF_000762265.1_ASM76226v1

Anyone know why this is happening?

I have experienced the same issue yesterday and after trying the --use-ftp argument without luck (this is proposed as an alternative in the forums), I have downloaded the database from Kraken 2, KrakenUniq and Bracken indexes (ref).


text:
How do I build the MinusB database for Kraken2? (Taxonomy issues)
I am attempting to build my own custom database for Kraken2. I have two questions:

If I have the MCPyV genome in a file called MCPyV.fasta, how do I build a database with just this?
How do I build the MinusB database? (https://benlangmead.github.io/aws-indexes/k2)

For number one, I attempted to build this by using:
kraken2-build --add-to-library MCPyV.fasta --db my_db

But was returned an error stating that there was no taxonomy information found. How do I make this MCPyV.fasta file amenable to Kraken2 while building the database?
For number 2, I again attempted to build the MinusB database using:
kraken2-build --build --threads 4 --db MinusB

and was returned an error regarding missing /taxonomy directory.
Does anyone know how to build these databases properly?

Hi @interestingQuestion61 the answer is helpful (I use Kraken too), please consider upvoting/accepting.

I am not sure if I understand your first question correctly.
Kraken is meant to be for the taxonomic classification of metagenomic sequencing data. Therefore, if you only have one target genome of interest there is no need to use Kraken. You can just use any sequence aligner and check the mapping rates.
However, if you are trying to add a custom genome to an existing database then you need a taxonomy id for the sequence names in the genome and add it to your DB using the --add-to-library command. If this is what you are asking this tutorial would help you https://bioinformaticsworkbook.org/dataAnalysis/Metagenomics/Kraken.html#gsc.tab=0
Regarding your second question, I truly believe those ref indexes are ready to go. You can download the tar file, extract it in a folder called MinusB and release the Kraken using, e.g., kraken2 --db MinusB yourreads.fq
Hope this helps


text:
How to Split VCF File into male and Female groups
Cross-posted on biostars
I have trio_example.vcf file and I want to split them into male and female groups. I tried this:
$ bcftools query -l trio_example.vcf

And it gave me:
N1
N2
N3

After that I tried this command:
for sample in `bcftools query -l trio_example.vcf.gz`; do vcf-subset --exclude-ref -c $sample trio_example.vcf.gz > ${sample}.vcf; done

to split 3 samples. Now its showing same data in all 3 N files
How can I solve this issue? Also I need to split the data into male and female. For splitting I am using:
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\n' trio_example.vcf | awk '{if($1 == "X" && $5 >= 0.3 && $5 <= 0.7) {print $6 " is male"} else if ($1 == "Y" && $5 == 1) {print $6 " is male"} else if ($1 == "X" && $5 == 1) {print $6 " is female"} else {print $6 " is of unknown sex"}}'

But it also not give me exact result

Do you know which of the individuals in the VCF are male and female? If so, you can use bcftools (view command) to subset the VCF twice, once creating a VCF only including males, a second time making a VCF only including females.

I dont know thats why I am using now this command: bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\n' trio_example.vcf | awk '{if($1 == "X" && $5 >= 0.3 && $5 <= 0.7) {print $6 " is male"} else if ($1 == "Y" && $5 == 1) {print $6 " is male"} else if ($1 == "X" && $5 == 1) {print $6 " is female"} else {print $6 " is of unknown sex"}}'

Is the information about the sex of your individuals anywhere in the VCF? In the INFO tag maybe? Somewhere in any paper or supplementary file? Which organism are you working on -- does this organism have sex chromosomes? Does the VCF have X and Y chromosomes (or Z/W etc depending on the sex determining system)? If so, individuals with variants at Y or hemizygous for X should be male, whereas individuals with homozygous variants at X should be female.

I only have X chromosome present in VCF File

What's the species. Drosophila?

I dont know about the species but when I separate N1 , N2 and N3 all files show the same record

If you don't even know what species you're working with I think you should figure that out ... because if you want to make an inference of sex from the VCF and that information simply isn't stated anywhere, you need to know the sex determining system you're working with. Drosophila only has X chromosomes, but males have one X and females have two. This kind of information is important if you want to infer sex because males can only have one variant at X, females can have two.

@AnasJamshed please add all extra information to your question, not in comments. As you can see your bcftools -f command is  very hard to read. Add it to the question please, along with all the other extra information you have given so that the next person to see the question can have all information there.

@arara don't forget the pseudoautosomal regions: you will likely have some, albeit few, homozygous variants in both X and Y even in males.

Cross-posted on biostars. OP was reprimanded on biostars by multiple users for opening a lot of tangential posts and is probably attempting a similar move here.


text:
Would duplicate climatic variables impact maxent?
I want to use maxent to make some future distribution predictions. I downloaded future climate data from worldclim. In worldclim's page there are three variables (max temp, min temp and precipitation) as well as a fourth titled "bioclimatic  variables" which mentions that it contains the other three and some other new ones. My question is: if i download all four and use them for my analysis, will maxent produce reliable results? Or would the duplicate variables be a problem?

When you say "maxent" and "distribution", what do you mean exactly? For example there is specific ecological software called "MaxEnt". But that's also a generic name for maximum entropy methods, which are fairly diverse.

Looking at the data source (which I believe to be this), it does not seem to line up with what you describe. There are more than 3 variables in each of the 4 datasets except the monthly values which has 3, but is presumably a bunch of time series. I am not sure that you will be able to meaningfully combine the time series data with e.g. bioclimatic point estimates for each location (?). You might think to yourself, "would it be handy to use a dataset that someone already preprocessed to be useful for biological modeling?"

Ok I understand better now so I'll write an answer.
Look at the bioclimatic data site. You can download them from the main page, but the suggestion to generate them from the time series data using an R package (dismo). The bioclimatic variables are thus various linear combinations of the 3 variable time series input.
Not having run this myself, I speculate, but it doesn't make sense to me that the bioclimatic variables are also time series, they appear to be summaries across the time series. So they're not exactly identical in information content, but they're close enough that I'm not sure it's worth the effort of figuring out how to munge the monthly data into a format that works with the summaries.
Someone already went to a lot of trouble to think about the biologically relevant representation of this data (e.g. what ANUCLIM does). I recommend relying on that work, or at the least looking at their arguments and deciding whether you agree enough to use their summaries.
It doesn't particularly matter what software you use, I'd argue (probably MaxEnt the species distribution package, but you haven't clarified). Most methods are going to be doing similar things with this data. They are going to find the variables with the best explanation and weigth them more heavily. Throwing in a couple extra correlated variables might get you a little margin, but it is unlikely to increase explanatory power.
If you want a formal methodology to compare the two, you could use something like AIC or BIC to compare the models formally with the 2 variable sets. But it hardly seems worth the effort of the data munging.
Assumption:
I am assuming that you are not talking about adding ~identical variables also into your model, e.g. files here, but are instead talking about the less identical monthly data here. I'm pretty sure this is right, but I wanted to clarify.