text:
enrichKegg arbitrarily excluding certain genes from gene set during analysis
I'm using the R-package enrichKEGG to perform an enrichment analysis on a vector of genes. As output I get a list of pathways containing some number of the genes in my vector. In the results there is a 'GeneRatio' column that shows how many of the genes in your vector appeared in each pathway.
My problem: My input vector contains 101 genes, but the GeneRatios for all the found pathways have 54 as the denominator (some of the ratios are: 12/54, 7/54, 8/54), suggesting that enrichKEGG is discarding nearly half of the genes in my vector before performing the analysis.
Ive tried adjusting every parameter: setting the pvaluecutoff and qvaluecutoff to 1, setting the minimum gene size to 1 and the maximum gene size to 10000, but its still 54 genes.
One possibility is that the ENTREZ-ids of some of the genes in my vector simply do not appear in any pathways. But I've looked at some of the genes that appear in the input-vector and don't appear in the list of 54 genes used in the analysis, and when I search those ids with (https://www.genome.jp/kegg/mapper/search.html) I find loads of results for pathways that the genes appear in.
Some of the genes that appear in the input vector AND in the list of 54 (ENTREZ ids):
"9319"   "10403"  "1026"   "9088"   "891"    "9700"   "890"    "113130" "9133"   "995"    "983"    "9212"   "3479"   "51512"

Some of the genes that appear in the input vector but DON'T appear in the list of 54 (ENTREZ ids):
"5980"   "3479"   "9319"   "51512"  "27092"  "27244"  "10403"  "8438"   "332"    "24137"  "4605"   "8411"   "2146"   "51655"

I'd really like to know based on what condition nearly half the genes are excluded from the analysis.

I'd suggest this is a ln2 cut-off which is not being cleared, this isn't dependent on probability but what is biologically reasonable before assessing a probability. This cut-off should be in the documentation. KEGG is very comprehensive db and 50% absence doesn't correlate with my experience, which even in massively diverse metagenomics data never rises above 40% against NCBI (i.e. worse case scenario). There's not really enough information to answer the question beyond this IMHO.


text:
Cellranger count error chnk
I run the cellranger count with this code:
cellranger count \
  --id=run_count_SRR11537950 \
  --transcriptome=refdata-gex-GRCh38-2024-A \
  --create-bam=true \
  --fastqs=SRR11537950 \
  --sample=SRR11537950 \
  --chemistry=fiveprime \
  --localcores=1 \
  --localmem=5

and error showed:
[error] Pipestance failed. Error log at: run_count_SRR11537950/SC_RNA_COUNTER_CS/SC_MULTI_CORE/MULTI_GEM_WELL_PROCESSOR/COUNT_GEM_WELL_PROCESSOR/_BASIC_SC_RNA_COUNTER/_MATRIX_COMPUTER/ALIGN_AND_COUNT/fork0/chnk0-u0de8b9d2e0/_errors
Log message: Job failed in stage code
signal: killed

What should I do? Thanks in advance!

I am almost certain that the localmem parameter is the issue here. CellRanger uses STAR aligner under the hood which is resource-hungry. We usually allocate 18 cores and at least 40GB of RAM to CellRanger while disabling all of its secondary analysis, such as clustering and tSNE. Five GB of RAM for sure are not sufficient.

I have to use upwards of 100GB to ensure things definitely run to completion.

@RamRS I just checked some older logs and we routinely use 18 cores and 40GB of RAM, but disabling all these secondary analysis CR can do (like tSNE etc). Still takes several hours to complete, but it usually works just fine.

We also use 40 GB RAM when running STAR with the human or mouse genome for regular RNA-Seq. It might be able to do with less memory if the parameters are tweaked, but that may not be possible when it's evoked through CellRanger.

source: Stack Exchange

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

sort: votes


text:
getGEO not working
trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE22nnn/GSE22093/suppl//GSE22093_RAW.tar?tool=geoquery'
Content type 'application/x-tar' length 357324800 bytes (340.8 MB)
downloaded 137.1 MB

Error in getGEOSuppFiles("GSE22093") :
  Failed to download E:/data/GSE22093/GSE22093_RAW.tar!

This error message appeared when I tried to download this dataset, I tried this function before and it worked, when I removed the files for some reason and tried again on the same PC and the same directory it didn't work,changed the directory didn't work, changed the permissions it didn't work out neither

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.

I could download the file from the above link ,its called GSE22093_RAW.tar. I suggest trying at a less busy time, e.g. weekend and considering whether there's a firewall issue blocking the download, e.g. download via tethering your mobile to you PC/laptop as a wifi hotspot. A firewall is the most likely explanation, NCBI do restrict/minimise bandwidth but not via ftp to my knowledge.

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


text:
Recent Bioinformatics Graduate Seeking Resources for NGS and Drug Discovery
I recently graduated with a bachelor's degree in bioinformatics. Unfortunately, the bioinformatics program at my university was limited, and I'm feeling the need to expand my knowledge in specific areas.
I'm particularly interested in next-generation sequencing (NGS) and  drug discovery.  I'd greatly appreciate any recommendations for resources that could help me gain a deeper understanding of these topics.
Specifically, I'm looking for:
Books: Introductory or intermediate-level texts that cover NGS technologies, data analysis, and applications in drug discovery.
Online Courses: Structured courses or tutorials that provide practical training in NGS analysis and its use in drug development pipelines.

I'm closing this question because it is of a general/discussion nature, and not related to a specific bioinformatics problem. Please ask discussion-type questions on discussion-based communities, like reddit.

@gringer I believe that reference request questions are not   general/discussion nature

For specific requests, yes, we have allowed it in the past (e.g. here). If you want to change to a specific topic, please do that, but I'm seeing 1) books on NGS technologies; 2) books on data analysis; 3) books on NGS use in drug discovery; 4) courses for NGS analysis; 5) courses for NGS use in drug discovery. It's better to ask multiple questions than to create a single post with multiple subquestions/areas. However, I expect you'll get a better response elsewhere with these requests.


text:
tabix errors when accessing 1000 Genomes data: "[E::bgzf_read] bgzf_read_block error -1 after 50219 of 52392 bytes" and Could not load .tbi/.csi index
I am trying to access 1000 Genomes (1KG) data using tabix as per the 1KG tabix documentation "How do I get a genomic region sub-section of your files?" but am getting errors.
Command I'm running:
I am getting errors when I run the sample command provided in the docs:
tabix -h ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz 2:39967768-39967768

Error message:
The error I'm getting is:
[E::bgzf_read] bgzf_read_block error -1 after 50219 of 52392 bytes
Could not load .tbi/.csi index of [...file name...]

Manual download
I tried to download a .tbi file and then run tabix as per this comment on GitHub but got an error
I ran:
wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/00-All.vcf.gz.tbi

which seemed to work, but this:
tabix ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/00-All.vcf.gz 1:17375-17375

returned an error:
hts_close returned non-zero status:

Other information
This worked for me a few days ago with the command above and my own queries, but I am now having issues.
I'm not very savvy with command line tools so it could very well be me, though the GUI Data Slicer tool  for accessing 1KG files doesn't seem to be working right now, so perhaps it's a server-side thing.
Software

I'm on a PC using the Windows Linux Subsystem
I installed tabix via bioconda after having issues compiling it directly as specified in docs
Since having this issue I have re-installed anaconda, bioconda, tabix, samtools, and htslib.

Tabix info
which tabix
/home/brouwern/anaconda3/bin/tabix

tabix --version

tabix (htslib) 1.3.1
Copyright (C) 2016 Genome Research Ltd.

try with https instead of ftp ?
$ tabix -h https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz 2:39967768-39967768  | cut -f1-8
##fileformat=VCFv4.0
##filedat=20101112
##datarelease=20100804
##samples=629
##description="Where BI calls are present, genotypes and alleles are from BI.  In there absence, UM genotypes are used.  If neither are available, no genotype information is present and the alleles are from the NCBI calls."
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth (only filtered reads used for calling)">
##FORMAT=<ID=GL,Number=3,Type=Float,Description="Log-scaled likelihoods for AA,AB,BB genotypes where A=ref and B=alt; not applicable if site is not biallelic">
##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GD,Number=1,Type=Float,Description="Genotype dosage.  Expected count of non-ref alleles [0,2]">
##FORMAT=<ID=OG,Number=1,Type=String,Description="Original Genotype input to Beagle">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=CB,Number=.,Type=String,Description="List of centres that called, UM (University of Michigan), BI (Broad Institute), BC (Boston College), NCBI">
##INFO=<ID=EUR_R2,Number=1,Type=Float,Description="R2 From Beagle based on European Samples">
##INFO=<ID=AFR_R2,Number=1,Type=Float,Description="R2 From Beagle based on AFRICAN Samples">
##INFO=<ID=ASN_R2,Number=1,Type=Float,Description="R2 From Beagle based on Asian Samples">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
2   39967768    rs11124691  T   A   .   PASS    DP=1991;AF=0.435;CB=UM,BI,BC,NCBI;EUR_R2=0.651;AFR_R2=0.742

Weird, not working for me same or similar error.

may be your proxy settings have changed ?


text:
Edge table on cytoscape won't show gene name on the shared name pp interaction column
Could you help me solve this problem?
So, basically both shared name and name column on this edge table not showing the gene name. It just show like this picture.
How can I change the name in the name table into gene name based on query term?

Please [edit] your question and clarify what you are expecting to see. Is this about seeing the names of the nodes in the network? Or about the table itself which is all you have shown us? Can you give us the actual data in text form (only a few lines, not the whole thing) so we can try to reproduce it? I don't see any gene names there, you seem to have protein names, I guess that's what you mean?

You need to use the mapSourceAndTarget tool (the little table icon in the header of the edge table) and map to display name.


text:
RNA virus (HIV, HCV, Influenza, or SARS-CoV-2) datasets with corresponding fitness values
Cross-posted on biostars
Does anyone know where I can find RNA virus (HIV, HCV, Influenza, or SARS-CoV-2) datasets with corresponding fitness values?
The dataset must contain genotypes and annotated .
I am interested in RNA viruses that affects humans .
Regarding to fitness it can be either of the following ( Epidemiological fitness, Adaptation Fitness, Immune Evasion Fitness, Transmission Fitness, Replication Fitness.
I want to visualise the dataset(s) using fitness landscape(s).
Thanks in advance!

Please [edit] your question and give us more detail about what you need. What would the dataset include? Genomes? Virus species names? Genes? And what do you mean by fitness? Replicative fitness? In what host? Transmission fitness? Are you only thinking of viruses that affect humans? Any RNA virus?

Thanks Terdon , for raising this , I have tried to give more info. I am still new to the field so my knowledge is limited.

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.

When you say "genotypes" what do you mean? And annotated how? With what information?

Genotype in this context I am refereeing to data representing the genetic makeup of the RNA viruses.
Annotations for the genotype provide additional context and information about the genetic data. These annotations can include functional impacts of mutations, related phenotypic data, references to known databases, and more.

I've worked on this stuff and a key part is identifying epistasis against selection pressure

You should sit down with a viral epidemiologist, immunologist, or virologist to identify the phenotypes and viruses you are targeting.
I personally would approach the question from phylogenetic "positive selection" analysis and then look at correlating phenotypes involving epidemiology, immunology and virology (in vitro replication). Again you need cross-functional support for this approach to be effective.
For SARS-CoV-2 the NCBI database contains all the relevant epidemiologic and virulence data in a singular (presumably SQL) database with a nice GUI frontend. It is really cool and thats the way to solve that virus. Check my post history for the exact location of the database on NCBI. With expert help with the interpretation of the vast outputs of NCBI's SARS-CoV-2 metadata thats absolutely doable in generic fitness landscapes, epistasis and all that. I have used NCBI SARS-CoV-2 database, not specifically in fitness landscapes/phenomena, but related genetic behaviour and it is fabulous. The speed this was unfurled in response to the epidemic was a credit to the NCBI.

Thank you very much for this detailed response

@user20708 thanks, however "accepts" and upvotes are much more appreciated than notes of thanks.


text:
LigPlot (LigPlus) Giving Different Interactions on Command Line and GUI
I am trying to determine the interaction of the targeted amino acid with other amino acids, but I am having trouble. When I pass a PDB file of a protein that has ligands, I can see the interactions. However, when I pass a PDB file of a protein without ligands, I don't see any interactions. Interestingly, when I run this ligand-free PDB file through LigPlot GUI, it shows the interactions.
PDB ID example: 1BBT
Here is the code I wrote:
filename = "1BBT.pdb"
ligPlotLoc = "LigPlus"
programLoc = f"{ligPlotLoc}/lib/exe_linux64/"
ligplot_prm = f"{ligPlotLoc}/lib/params/ligplot.prm"
atom_number = 82
chain_ID = 1

os.system(f"{programLoc}hbadd {filename} {ligPlotLoc}/components.cif -wkdir ./")
os.system(f"{programLoc}hbplus -L -f hbplus.rc -h 2.90 -d 3.90 -N {filename} -wkdir ./")
os.system(f"{programLoc}hbplus -L -f hbplus.rc -h 2.70 -d 3.35 {filename} -wkdir ./")
os.system(f"{programLoc}ligplot {filename} {atom_number} {atom_number} {chain_ID} -wkdir ./ -prm {ligplot_prm} -ctype 1")

Below is the image of the interaction I obtained from running the LigPlot command line.

Below is the interaction that I am trying to get using the LigPlot command line but haven't been able to. However, I successfully obtained it using the LigPlot GUI. But wanted to get it done through ligplot command

its not free, so you better ask to the user community or here https://www.ebi.ac.uk/thornton-srv/software/LigPlus/contactus.html In case of any problems or questions relating to LigPlot+, please contact (preferably by e-mail):
Roman Laskowski
European Bioinformatics Institute,
Wellcome Trust Genome Campus,
Hinxton,
Cambridge, CB10 1SD,
United Kingdom
Phone Number: +44 (0)1223 492 542
Fax Number: +44 (0)1223 494 468
Email: roman@ebi.ac.uk
Personal Page: http://www.ebi.ac.uk/~roman/

Thank you @pippo1980


text:
Can't convert my vcf to plink
I have run this and checked the .map for number of contigs and specified in the code plink --file vcf_corrected --chr-set 897 --make-bed --out output_file but keep getting the error below
PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to file.log.
Options in effect:
  --chr-set 897
  --file vcf_corrected
  --make-bed
  --out outputfile

Error: Invalid --chr-set parameter '897'.
For more information, try "plink --help <flag name>" or "plink --help | more".

Please provide more context around your question. Are you trying to use Plink on a draft genome assembly with 897 contigs?

@gringer: Any idea how my comment disappeared? Asking because you're a mod

@RamRS answers should be written as answers, not comments. See here.

Update
Thanks to @RamRS for pointing out: the documentation that I quoted below is apparently a sort of inside joke (!) for the plink devs.
It sounds like you might be able to use the flag --allow-extra-chr to support a draft genome:

Contact us, and we'll send you a build supporting a higher autosome limit. Note that this isn't necessary if you're dealing with a draft assembly with lots of contigs, rather than actual autosomes—the standard build can handle that if you name your contigs 'contig1', 'contig2', etc. and use the --allow-extra-chr flag."

Original answer
As @gringer suggests, I think that the plink docs tell you what you need to know:

--chr-set changes the chromosome set. The first parameter specifies the number of diploid autosome pairs if positive, or haploid chromosomes if negative. (Polyploid and aneuploid data are not supported, and there is currently no special handling of sex or mitochondrial chromosomes in all-haploid chromosome sets.)
Given diploid autosomes, the remaining modifiers let you indicate the absence of specific non-autosomal chromosomes, as an extra sanity check on the input data. Note that, when there are n autosome pairs, the X chromosome is assigned numeric code n+1, Y is n+2, XY (pseudo-autosomal region of X) is n+3, and MT (mitochondria) is n+4.
n is currently limited to 95, so if you're working with adder's-tongue fern genomes, you're out of luck.

I don't think that plink will work with your current assembly.

I had added a comment but it looks like it's disappeared somehow. The docs specify the "you're out of luck" as a joke. Read further: "Just kidding. Contact us, and we'll send you a build supporting a higher autosome limit. Note that this isn't necessary if you're dealing with a draft assembly with lots of contigs, rather than actual autosomes—the standard build can handle that if you name your contigs 'contig1', 'contig2', etc. and use the --allow-extra-chr flag."

Thanks @RamRS. Apparently I should take my own advice on reading the docs! However I'm a little surprised that even as a joke they are misinforming the reader.


text:
How to only enrich sensible GO and KEGG terms?
I'm thinking of analyzing enriched pathways using KEGG and GO databases on de-novo assembled transcripts using clusterProfiler 4.0. After annotating them with eggnog-mapper, I have some annotated and some unannotated Transcripts. Here is the doubt:
I'm only interested in pathways and ontologies that make sense (Like "Lipid Biosynthesis" or "Antibiotic degradation", not like "Metabolic process" or "Cytosol"). So I was thinking, after getting the p-values (not adjusted yet), I will only do post-hoc analysis for gene groups I care about.

Is this strategy sound?
Can it be automated somehow? because the list of ontologies and KEGG modules is so long? Because I heard somewhere that there is a parent-child graph-like relationship among GO terms. Also, a recent update in clusterProfiler 4.0 performs this analysis on KEGG modules, but I don't know how to deal with the useless portion of GO terms.

Thanks for your kind help.


text:
Facing problem with BLASTx using BioPython in Google Colab
I'm trying to BLASTX a sequence using Biopython in Google Colab. Here is the code I am using. I'm not able to understand why I am getting an empty list [] for alignments.
result_handle = NCBIWWW.qblast("blastx", "nr", test_seq.seq, expect=1e-05)
blast_record = NCBIXML.read(result_handle)
alignments = []
for alignment in tqdm(blast_record.alignments):
  for hsp in alignment.hsps:
      alignments.append({
          "title": alignment.title,
          "length": alignment.length,
          "e_value": hsp.expect,
          "score": hsp.score,
          "query": hsp.query,
          "match": hsp.match,
          "subject": hsp.sbjct
      })
len(alignments) # its coming 0 because its empty

As a note, here is the test sequence I am using, which is giving proper results in BLASTX on NCBI server in NCBI NR database. Any help would be very helpful.
>test_seq
GCAGCAGCAGCAGCAGCAGCGCGGCGCGCGACTGCCTGGCCATGAATGCTTTCTCCAAGCTGCTGCGCGT
GCTGACCACGCCCTCCATGCTGCACAACCTGTGCAACAGCGACCAAGGCGTGCACGTGCTGATCCGGCTC
ATCAAGCTGAGCTCGCCGCCAGAGACGGCTTTCATTGCTGACCAGGTCCTGCACGACTTGGCTGACATGG
TGGTCAGCAAGGGCGGCGCTCTGGTGGTCTGCACCCTGGTCGACCGCCTGGCCGTCGCGGCCAGCTGCGA
CAACCCGGGGGTGGTGGAGCGGGCCACGGGCGCGGTTGTGCGCATGTGCAGCCACCTGTGCCAGCACCAT
CGCGCGCTGCAGCAGTGTGCGCGGCACCAGTACCTGGGTCCCAAGGTCGTTGAGCTGCTGCAGTGCGGCC
TGCCCGCTGTGGATGCCTTCAAGATGGGGTTTGAGCTGGTCCAGAACGCGCCCGCGCTGGCGGCCACCCG
CCACGGCAGCCGCACGCTGCAGCTGCTCACGGAGATTCGGGCCGCCGATGCGCGGTCGCGAGAGATGCTG
CGAGACGTGGTCTGCCACGCGTCCAAGGCGCTGCAGGGCCACTTTGGCCAGCTGGCGGCCAAGGGGGAGA
CGGGAGGGCCCGCCCTGGTCCGCGCCCTGCTGGACTACCTGCTGGCGGAGCGGGAGGACGACTGGGTGGT
GCAGATCACCTGGGAGCTGCTGGCGGCGGGCAGTGAGCTGCGGCGGGACCGGGAGGCACATGACTTGCTG
GTGCACAGCCTGGCCCTGCCGCTGTGGCGCCCCACCGAGGCGCAGCAGCACTGGGCGGCGCTTGCGGCCG
CTGGCGTGGATGTGGCAGCCCTCAAGCCCAAGGTGGATGCAGCCAAGCAGCACCAGCCGCGGCTTCCGGA
GCGGGTGCTGGCATCGCTGTCCGAGGTGCAGCCTTCGTCACCTGCAGAGGTCAGCGCCAAGCAGCGCACG
AGCTGTCTTCCCCGCTCCGCATGCTCAGGGCCGGGCTGCTCCTGCCTGCCCGGAGTGGTGAGCGGGATCC
TGACGTGCCTCCCTGCCCCTCGCTTGTCCAATGATGTTCTATATTCCATCGTCGCAGGCTGTGCGGCGCT
TCCCTTACTATCTGCGGGATGGAAGCCAGGCGGCTCGCAACTCTGTGCTCAGAAAGCAGCAGGGCGCAGG
GCAGCCAGCCGTGCCCCCGCCGCCCCCGCCGCGGGAGCAGCAGCAGCCCAAGGCGCCCCCTGGCTTTGGG
GCGGCAGGAGCCGCTGCAGCGCCAGCGGCTCCTGCGGCTGCCAAGGTGCCCCCACCCGGCTTTGCGGCAG
CGCCGTCGGCGGCGCCAGCAGCAGCGCCAGCGGCCATCACTTTCCCAAATGGCGCAAGCGTGCCCATGCC
GCAGCCCTATGGCGCCAGCAGCAGCAGCGCTGCGCCTGCAGCCAGCAAGGCGGTTCCGCCGCCAGCCCCA
CAGCAGCCAGCTGTCGTAGCCGCTGCTGCGATGGTGCCCCCGCCCTCGCCCCTCAAGCCTGCAGCCCCTC
CGCAGCAGCCAGTGGCACTGCCACCGCCTGCCCCCATGCAGCCAGCAGCAGCGCCCCTGGTGCCTGCACC
TGTGCCGCAGCCCGCGCCTGTGCAGCCACCAGTGGCAGCGCCCCGGCCCCCGGCACCTGTGCCGCAGCCT
GCTGTGCAGCCACCTGTGGCGGTGCGGCCACCAGTGCCACCACCTGCGCCTGTGGCACCCGGCATGAATG
CAGCAGCCCGGCCATTTGTGCAGGCAGTGCCTGCCCCTGCGGCGCAACCCCTGCCCCCGCACCTGCAGCC
CAACCTGCAGCAGCAGCAACAGCCGCAGCAGCCCGCGCACCTTCAAGCGCTGCCGCCGCACCTGCGGCAG
CCGGTGGCAGCGGCAGCAGGCAGCAGTGCCGGGCAGGCGGGCACGGCAGCACCCATCACGCCTGTGCCGG
CACAGCACCCTGCGGGCGGCTTCACCGCGCCGCCTGCACCGCAGCAGCTGCTGCAGCAGCTGCTTGGCAG
CGCTGCCCCGCGGCCCGCGGTCCCGGCAGCGGCGCCAGCAGCGGTGCGGCCGCTTGCTCCGCCGCAGGCA
GCAGCAGCTCAGGGCAGCGCAGCTGCCTTTGATGGCTACGGCAGCCTCACAGAAATGCTGTCTGAGCTAA
CTGCCGCTGCGGCGCCACAGCCGGTGGCACCTGCACGCCCAGCCTACCGCCCGGCCTACCGCCCGCCTGG
AGCAGCTGCTTACCGCCCGCCCACCGCCGCGGCGCCGCCGCCACCGCCCGCGCCAGCGGCCGCACCTGCG
CCTGTTTCTGTGCCTGTGGGCATTCCAGCCCCCTCGCCTGCAACCCGGTTCCCCATGCCCGGCGCGCCAC
GGCCCGCACCGGTGCCAGCAGCTGCGGCGCCGCCGCCGCCGGCAGCGCCCGCCCGCCCGGCGGTGCCAGC
GGTGCCGGTGGCATCGGCGACCAACTGGACCTGCAGGGTGTGCACCTACCGGCATGAGGGAGCGGAGGCT
CAGTTCCTGTCCTGTGCCATGTGCGGCTCCCCACAGTCCTAGGCAGCCTGCTGTGGCAGTCCCTAGGCTG
GCTGCTGTGGCAGTCCTGGGCAGGCCGCTTGGTGGCTAGCTGGCTTCGCTGGCTCGGCAGCACTTCAACG
CTCCCACAACCTTCCATTGTTCCTTGCACTGCAACCCTTGCCCCCTGATCGCATCGCTCGCCCTTGCCTC
TTGTCATTCATGCTGCACCAACCCAGCCCCTGTGCATGTAACACCAGAGATGACCCCCCTGTTGTATTGG
GCTCGCTTGCTTTGAGGCACGGGGGTTGTTTCACAGGAGGAATGAACAGGGGGGCACATGGGACGGCTGG
GAAGAAGCAGCACAAAGAGCAACAGGAACAGAATCATCAAATGAGCATGTGGTGGTACACTTGACTTCGG
GGGATGAAGAAGAAGAGGCAGCCAGGGGTGCTGGGGCGGAGGGGGGTGTCGATCCGAGGCATCATGCCGC
CAGCGGTAAGGGTGCAAGCAGTGGGGGTGCGTATTTAGACAGCAGGGCATTCGTGGGGCAAGGAAGAGAT
GGAGGAGTGTGGCGGGGCTGGCAGGGGGTTCAATGCAGCCTTCCTGCCGCACACTTTGCAGACATGCATT
TTCCAGGGAAGAGACGAGAGGAAGAGCAGTGCGACGGGGCGGGACAGCGAAGGGGCGGGCGCACGGGGGC
GCAGTCCCGCCCATGGCGCAGCCGCCACGGCCGGAACAAACATTTTGTTTTTGATGGTACACAAACAAAC
AAACAGGAGCGATGTTAAGAGTTAGGGAGCGGAGCAATGCAGAGGTGGGTGGGCAGGTAGCAAGGCAGGG
AGCGGGTATATGCCGGGGGCGGGCCACAGTGCAGTGGCGGCATGGCGGCAGCGCCACGCTGCCGCCGCCC
ACACCGACGGCTGGGCGCCGCCACTTCAGTCGGCGAGCGTCGCGCGCGCCAGCTGCGCTTGCACAGCCGC
CGTGCCGCCGCCCTGCTCGCCAGATGCAGTGCGCACCGGCTGGGCAGCGGAAGCAGACTCCTCAGACTGC
TCCAGAGAGCTGGCAGAGTCCACCTTGCCCTTGCTGCGGCTGCCGTTGCGGCTCAGCATGTGCAGCGGCC
TGCGGCGAAGCAAGGGCGGTGCGGCGGGGGCCGCAGCCAGGGTCGGCATCCAAGCCGCAGCCGGCAGCAG
CGTGCGGTACAGAGGGCGGCTGGCCTCGGCCTGCCGGCGCTGCCAAGCATGGGTGCTATGAGGAAGCAGG
CCTGCAGAGTTGCAGGCCGGCAGGCAGGGCAGCGCCAATCAGGGCGGGCACCAGCAGTGCCAGGCTCGCC
AGCAACAAGCTGCAGCCCTTTGCCACACTCACCGCGTGTACGGCGGCTGCGACTGCTGCGACTGCTGCTG
CGCAAACGGCGTGCCGCTGTAGTGGCCGCCGCCGCCGCCGCCAGCGGGGTGCTGCGGCGGGAATGCCGGC
GATGCCATGCCGGTGCTTGTGGGGTAGCCGTAGGCAGCTGCGCCGGGCACATAGCCGCCCATCTGGCCTC
CCATCATCATGTACGGGCTGGGGGCCATGCCCCCCATGCCAGCCATTGGGCTCATGCCACCGTACTGCTG
GGCCATGCCGGCATACTGCATGCCAAAGCCGGGGCCGCCCGGCACCATGCCTGCATGCGGGTAGCCGAAC
ATGCCAAAGCCGGCGGCCGCGTAGGGGCTGGGCGGGCTGCCAAACAGCGCGCCGGGTGAGCCCACCACGG
CGCCGGGCGGGCCGCCAAAGGGCGGCAGGCCCATGCCGCCGCGGCCGCCGGCGCCGCTGGTGCTG

EDIT: By the way, to test the code, I tried running Biopython in Jupyter locally. When I tried with BLASTp with same code for a different sequence, I got this error:
URLError: <urlopen error [Errno 11001] getaddrinfo failed>

I'm thinking there is some problem with the scraper implemented in Biopython. When attempting BLASTX, the code keeps on running for an unreasonable amount of time (>10 min. I had to stop the Kernel). Let me know if it makes sense.
Thanks

Can you show us the complete program, including any include directives so we can reproduce this locally? (not sure why you chose to undo my edit, by the way, that is really hard to read)

Is that really your sequence? You have everything, including the header on a single line? And should we assume test_seq.seq is a file? What kind of file? Simple text file? Windows line endings? Unix?

test_seq.seq is definitely a Bio.Seq.Seq object. I checked that 5 times before posting this query. NCBI servers seem to be having a problem in general, really.

Regarding the code... it's just Biopython all the way. From FASTA parsing to getting the first sequence.

That's not particularly helpful. If you're not providing at least stub code for others to see the same problem, they're less likely to help you out with a solution.

I've been able to replicate the issue you've had (or at least produce a similar issue) with the following code:
#!/usr/bin/env python

## initial setup code from 'www_blast' biopython example
## (Brad Chapman / Biopython)

# standard library
from io import StringIO

# biopython
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

from tqdm import tqdm # progress bar

# first get the sequence we want to parse from a FASTA file
test_seq = next(SeqIO.parse("test_seq.fasta", "fasta"))

# this fetch query takes a *long* time
result_handle = NCBIWWW.qblast("blastn", "nr",
                               test_seq.format("fasta"),
                               expect=1e-05)

# save the results for later, in case we want to look at it
with open("m_cold_blast.out", "w") as save_file:
    blast_results = result_handle.read()
    save_file.write(blast_results)

blast_record = NCBIXML.read(result_handle)

alignments = []
for alignment in tqdm(blast_record.alignments):
  for hsp in alignment.hsps:
      alignments.append({
          "title": alignment.title,
          "length": alignment.length,
          "e_value": hsp.expect,
          "score": hsp.score,
          "query": hsp.query,
          "match": hsp.match,
          "subject": hsp.sbjct
      })
len(alignments) # result: 0

Looking at the message in the saved XML output - I'd recommend doing similar with saving and looking at intermediate data when debugging issues - I see the following:
<Iteration_message>[blastsrv4.REAL]: Error: CPU usage limit was exceeded, resulting in SIGXCPU (24).</Iteration_message>

So my guess is that your query string is too long.
After further investigation, I've come up with the following code, which is able to properly process a shorter string (280bp):
#!/usr/bin/env python

## initial setup code from 'www_blast' biopython example
## (Brad Chapman / Biopython)

# standard library
from io import StringIO

# biopython
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

from tqdm import tqdm # progress bar

# first get the sequence we want to parse from a FASTA file
test_seq = next(SeqIO.parse("test_seq_short.fasta", "fasta"))

# this takes a long time
result_handle = NCBIWWW.qblast("blastx", "nr",
                               test_seq.format("fasta"),
                               expect=1e-05)

# save the results for later, in case we want to look at it
with open("m_cold_blast.out", "w") as save_file:
    blast_results = result_handle.read()
    save_file.write(blast_results)

string_result_handle = StringIO(blast_results)
blast_record = NCBIXML.read(string_result_handle)

alignments = []
for alignment in tqdm(blast_record.alignments):
  for hsp in alignment.hsps:
      alignments.append({
          "title": alignment.title,
          "length": alignment.length,
          "e_value": hsp.expect,
          "score": hsp.score,
          "query": hsp.query,
          "match": hsp.match,
          "subject": hsp.sbjct
      })

>>> len(alignments)
4
>>> alignments
[{'title': 'gb|PRW60132.1| 60S ribosomal L15 [Chlorella sorokiniana]', 'length': 1087, 'e_value': 2.9921e-51, 'score': 466.0, 'query': 'SSSSSSAARDCLAMNAFSKLLRVLTTPSMLHNLCNSDQGVHVLIRLIKLSSPPETAFIADQVLHDLADMVVSKGGALVVCTLVDRLAVAASC', 'match': 'SSSSSSAARDCLAMNAFSKLLRVLTTPSMLHNLCNSDQGVHVLIRLIKLSSPPETAFIADQVLHDLADMVVSKGGALVVCTLVDRLAVAASC', 'subject': 'SSSSSSAARDCLAMNAFSKLLRVLTTPSMLHNLCNSDQGVHVLIRLIKLSSPPETAFIADQVLHDLADMVVSKGGALVVCTLVDRLAVAASC'}, {'title': 'gb|KAI7846133.1| hypothetical protein COHA_000394 [Chlorella ohadii]', 'length': 864, 'e_value': 1.83868e-41, 'score': 394.0, 'query': 'AARDCLAMNAFSKLLRVLTTPSMLHNLCNSDQGVHVLIRLIKLSSPPETAFIADQVLHDLADMVVSKGGALVVCTLVDRLAVAASC', 'match': 'A  DCLAM +FSKLLRVLT+ SMLHNLC SDQGVHVL+RLIKLSSPPETAFIADQV+HDLADMV SKGGALVVCTLVDRLAVAASC', 'subject': 'AGHDCLAMTSFSKLLRVLTSASMLHNLCMSDQGVHVLVRLIKLSSPPETAFIADQVMHDLADMVTSKGGALVVCTLVDRLAVAASC'}, {'title': 'ref|XP_005845358.1| hypothetical protein CHLNCDRAFT_137175 [Chlorella variabilis] >gb|EFN53256.1| hypothetical protein CHLNCDRAFT_137175 [Chlorella variabilis]', 'length': 1059, 'e_value': 2.1613e-11, 'score': 171.0, 'query': 'CLAMNAFSKLLRVLTTPSMLHNLCNSDQGVHVLIRLIKLSSPPETAFIADQVLHDLADMVVSKGGALVVCTLVDRLAVAAS', 'match': 'C A+ A+S+LL+ LT P+M H L  SD+GV+VL++L+KL SP ET +IA+Q    + +   SK G+L V  LVDRLA  A+', 'subject': 'CTALRAYSELLKALT-PAM-HQLATSDRGVYVLVKLVKLGSPMETVYIAEQATPMVVEWSGSKDGSLFVTQLVDRLAQVAT'}, {'title': 'gb|KAI3436075.1| hypothetical protein D9Q98_002133 [Chlorella vulgaris]', 'length': 1091, 'e_value': 4.23077e-06, 'score': 132.0, 'query': 'CLAMNAFSKLLRVLTTPSMLHNLCNSDQGVHVLIRLIKLSSPPETAFIADQVLHDLADMVVSKGGALVVCTLVDRLAVAAS', 'match': 'C AM A+S+LL+ L     +HNL   D+GVH L++L+   S  E  ++A+Q    + D + SK G L    L+ RLA +A+', 'subject': 'CTAMMAYSELLKALAL--AMHNLSMRDKGVHALVKLVDTGSTLEAVYLAEQATPLMLDWMGSKDGILFAVKLIHRLAQSAA'}]

Comparing with the code you have provided, the only difference I can see is that my working result text is wrapped inside a StringIO Structure.

I tried BLASTp using a 46aa sequence. That didn't work either. When I tried a 25-sequence FASTA file in NCBI, there was the same error (CPU usage exceeded) I guess, NCBI servers are getting overwhelmed...

Yeah, I can't imagine that the query is too long. That isn't a particularly large sequence, so if it couldn't handle that, the whole blast tool would be useless.


text:
Should I include unannotated transcripts in a pathway enrichment analysis of a de-novo assembled transcriptome?
I'm trying to perform an enrichment analysis on a set of Transcripts assembled de-novo, using a GSEA like method on KEGG and GO. The set contains annotated and unannotated transcripts. My question is:

Should I include the Transcripts which don't have a GO and/or KEGG term?
Will it matter whether it is included or excluded?

Thanks for your kind help.


text:
pca plot from multiple sequence alignment
I would like to get a pca plot from the multiple sequence alignment (fasta file). I extracted snps in vcf format using snp-sites program.
snp-sites -v -o output.vcf alignment.fasta

SNPRelate package was used to get a pca plot from the vcf file
library(gdsfmt)
library(SNPRelate)

vcf.fn <- "C:/users/downloads/output.vcf"
snpgdsVCF2GDS(vcf.fn, "test.gds", method="biallelic.only")
snpgdsSummary("test.gds")

The total number of samples: 266
The total number of SNPs: 0
SNP genotypes are stored in SNP-major mode (Sample X SNP)

genofile <- openfn.gds("test.gds")
test_pca<-snpgdsPCA(genofile)
Hint: it is suggested to call `snpgdsOpen' to open a SNP GDS file instead of `openfn.gds'.
Principal Component Analysis (PCA) on genotypes:
Excluding 0 SNP on non-autosomes
Error in .InitFile2(cmd = "Principal Component Analysis (PCA) on genotypes:",  :
  There is no SNP!

How do I resolve this error? or Any other methods available to get a PCA plot from the alignment file?

What is the error? You have no SNPs in your input so cannot do any analysis on SNPs. For more help, you will need to edit your question and explain what you are doing in more detail. What are you aligning here? Do you have any data in the vcf?

Please show the first few lines of your VCF file, especially the tiles that include genotypes.


text:
Unable to install Enrichment app
I am unable to install the enrichment app in my cytoscape 3.10.2.
Undefined error was shown on the screen. I have downloaded the Enrichment map app 3.4.0 version.
Can someone please suggest me which jar file I have to use to install the app.
Thank you in advance

Please [edit] (don't answer in the comments) your question and i) tell us your operating system, ii) what this "app" is, where you got it, and exactly how you tried to install it so we can try to reproduce what you saw.

Do you mean the EnrichmentMap App? https://apps.cytoscape.org/apps/enrichmentmap

Please copy-paste the text of the error, or (if that is not possible) a screenshot of the error.

Please provide more information.  EnrichmentMap seems to install fine for me, but I'm on a Linux machine.  What type of OS are you using? What other apps do you have installed? How are you attempting to install it?  All these would be helpful.


text:
How to interpret mutagenesis results from pymol to make conclusions about structural/functional changes and predict resulting pathogenesis?
I'm a beginner bioinformatics student. I am currently working on my first project. We have gathered data about SNPs that are associated with a particular disease. We aim to understand the functional and structural implications of these SNPs. I have run mutagenesis on Pymol and gathered data on the variation in Hydrogen bonds.
We are now unsure how to analyze this data or interpret the mutagenesis result to get to a point where we can say that:

A particular structural change is occurring and
This change is now leading to functional changes thereby making it pathogenic.

Is it just basic understanding? Or are there any other software tools to help with this interpretation?
Thanks in advance!

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.

Is this human data?

Assuming this is human data, this seems like a job for the Ensembl Variant Effect Predictor (VEP):

VEP determines the effect of your variants (SNPs, insertions, deletions, CNVs or structural variants) on genes, transcripts, and protein sequence, as well as regulatory regions.
Simply input the coordinates of your variants and the nucleotide changes to find out the:

Genes and Transcripts affected by the variants
Location of the variants (e.g. upstream of a transcript, in coding sequence, in non-coding RNA, in regulatory regions)
Consequence of your variants on the protein sequence (e.g. stop gained, missense, stop lost, frameshift), see variant consequences
Known variants that match yours, and associated minor allele frequencies from the 1000 Genomes Project
SIFT and PolyPhen-2 scores for changes to protein sequence
... And more! See data types, versions.

If you can get your data into VCF format, you should be able to run VEP on it, which will tell you if there is a protein change, and whether that change is predicted to have a biological consequence.
Aside on terminology: in bioinformatics, when the term structural change is used, it usually refers to large-scale genetic changes (e.g. large deletions, translocations, insertions, or inversions). As SNPs are single point variations, they don't pass the threshold for structural change.

Thank you so much !


text:
Protein conservation analysis with trident scores in R script
I am currently working on analyzing protein conservation across amino acid positions using Trident scores as output from the MstatX tool. I have 40,000 Trident score files, each containing amino acid positions and their respective scores.
To streamline my analysis, I initially wrote a bash script to loop through each file and then implemented an R script for the calculations. The R script calculates moving averages with different spans to identify regions with high deviation in local protein conservation, referred to as "hotspots." These hotspots are then analyzed to determine the strongest deviation, its position, and related metrics.
Here is a snippet of my R script:
# assign input files to a vector
args = commandArgs(trailingOnly = T)

# put sample name in variable
sample<-sprintf("%s", args[1])

# load sample and make data frame
raw.entropy<-read.table(sprintf("%s", args[1]))
col.names<-c("AA.pos", "Entropy.Score")
colnames(raw.entropy)<-col.names
AA.pos<-raw.entropy$AA.pos
Entropy.Score<-raw.entropy$Entropy.Score

# make sample name nicer
sub("_entropyscores", "", sample)->sample

# load ggplot
library(ggplot2)

# calculate short and long moving averages
mydata<-ggplot(data=raw.entropy, aes(x=AA.pos)) +
  geom_point(aes(y=Entropy.Score), size=2, color="black") +
  geom_smooth(mapping=aes(y=Entropy.Score), method = "loess", color = "red", span = 0.1) +
  geom_smooth(mapping=aes(y=Entropy.Score), method = "loess", color = "blue", span = 0.4) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = sprintf("%s", args[1]), x = "Amino Acid Positions", y = "Trident Score")

# infer ratios between short and long span moving averages
# ratios represent extent to which local conservation is adhered to
# ratio around 1 means perfect conservation, ratio around 0 means poor conservation
rolling_mean_10<-layer_data(mydata, 2)$y
rolling_mean_50<-layer_data(mydata, 3)$y
ratios<-rolling_mean_10/rolling_mean_50

# find index value of 10th percentile
lowerbound<-NULL
for (i in (0.1*length(rolling_mean_10)):length(rolling_mean_10)){
if (rolling_mean_10[i] >= 0.9*median(raw.entropy$Entropy.Score)){
lowerbound<-i
break
}
}

# find index value of 90th percentile
upperbound<-NULL
for (i in (0.9*length(rolling_mean_10)):1){
if (rolling_mean_10[i] >= 0.9*median(raw.entropy$Entropy.Score)){
upperbound<-i
break
}
}

# find location in sequence with highest deviation in local protein conservation
which(ratios==min(ratios[lowerbound:upperbound]))->hotspot
ratios[hotspot]->cutoff
layer_data(mydata, 3)$x[hotspot]->x
percentile<-x/length(raw.entropy$AA.pos)*100

# write relevant items to standard output
# order is cluster name, strongest deviation, rounded amino acid position, percentile where deviation occurred
cat(sprintf("%s %s %s %s %s\n", sample, cutoff, round(x), round(percentile, 2)), file="./output_deviation.txt", append=TRUE)

After writing the output file, I aim to create a density plot with the x-axis representing the strongest deviation of each input file. However, the current method does not capture the expected bimodal distribution. I am seeking advice on alternative ways to calculate small but strong deviations beyond using simple moving averages.
Some options I've explored include EMA, z-scores, step detection, and anomaly detection. I would greatly appreciate any insights or suggestions.


text:
Comparing sc-seq combination and bulk-sequencing
I have expression matrices for different cell types, representing the expression of individual cells of that type. They were learned through a generative model, so I am confident they represent the approximate expression patterns of specific cell types. Now, I want to implement a bulk sequencing deconvolution using the aforementioned expression patterns. The pseudo-bulk I used is the summation of a large number of single-cell sequencing data (sc-seq). That's the background.
My first approach is to design an optimization process to optimize a series of weights, so that the product of the weights and the cell type expression approximates the pseudo-bulk. I was advised to use Poisson loss as the loss function because it aligns with the biological characteristics of RNA-seq. However, I couldn't add non-negativity constraints to the weights during optimization, resulting in negative values in the optimization results, which is meaningless. Then I found the optimiza.nnls method in the scipy package, which implements non-negativity constraints, but it uses Euclidean distance to compare the pseudo-bulk and the sc-seq combination. I obtained some good results using this method, but I have the following questions.

Can I use Euclidean distance to compare the differences between two sequencing methods? To me, this problem seems to become a linear regression problem, i.e., combining sc-seq to approximate pseudo-bulk. At this stage, there don't seem to be any biological distribution assumptions, so I guess it's feasible.
If the answer to the previous question is no, what biological assumptions does using Poisson loss follow, and what am I ignoring when using Euclidean distance for comparison?
If I want to continue using Poisson loss to optimize weights, how should I set the non-negativity constraints on the weights? I have tried methods such as ReLU and softmax in machine learning, but the results are not good.

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

I do not know the specific answers to your questions but I know there are several published papers out there about single cell deconvolution that might help.


text:
How to import data from cell ranger to R (Seurat) using Read10X?
I am loading cellranger data into seuratobject as following, but it keeps showing Error in
Barcode file missing. Expecting barcodes.tsv.gz

the running code and result:
#clean up your environment
rm(list=ls())

# Load required libraries
library(Seurat)
library(dplyr)

# Load required library for reading Matrix Market files
library(Matrix)

# For output from CellRanger >= 3.0 with multiple data types
data_dir <- 'C:/Users/Administrator/Desktop/10x_raw'

list.files(data_dir) # Should show barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz
[1] "GSM7494257_AML16_DX_raw_barcodes.tsv.gz" "GSM7494257_AML16_DX_raw_genes.tsv.gz"
[3] "GSM7494257_AML16_DX_raw_matrix.mtx.gz"

data <- Read10X(data.dir = data_dir)
Error in Read10X(data.dir = data_dir) :
  Barcode file missing. Expecting barcodes.tsv.gz

seurat_object = CreateSeuratObject(counts = data$`Gene Expression`)
Error in data$`Gene Expression` :
  object of type 'closure' is not subsettable

I am sure this file inside the directory. Dose anyone know what's the issue?
Thanks in advance.

The list.files line demonstrates that the precise files it's looking for are not there. What do you mean by "I am sure this file inside the directory"? Can you show a listing of files in the directory that you are trying to read from?

It is been a while since the last time I have done this, but the problem seems to be your filenames. For example Read10X() is expecting a file barcodes.tsv.gz but you seem to have sample prefixes in your file names GSM7494257_AML16_DX_raw_barcodes.tsv.gz.
If I remember correctly, cellranger outputs a directory sth like sample-name/outs/ and this outs dir contains barcodes.tsv.gz, genes.tsv.gz, matrix.mtx.gz that Seurat can easily read.

Thank you much, you are right! Exactly the name without prefixes.


text:
How to visualize unique and overlapping genes in upSet R plot
I am still a newbie in the field of bioinformatics. I have been attempting to visualize unique and common genes across six different tissues using the UpSetR plot package. This script was used to generate the plot; however, out of the six tissues, only six could be on the plot.
library(data.table)
library(UpSetR)
d <- data.table::fread(text = "Brainstem    Cerebellum  Olfactory Bulb  Heart   Kidney  Liver
Rac3    Coq9    Ndufs7  Slc25a10    Aacs    Acss3
Ndufb5  Mtnd5   Cotl1   Adhfe1  Aass    Cyp4v2
Clic6   Kpna4   Gabarapl2   Ddx5    Abcc4   Slc7a2
Rbmx    Snrpg   Rbm3    Pgd Acot12  Ehhadh
Ptcd3   Micos13 Atp5f1e Gss Acsf2   Selenbp2
Ndufa13 Ndufa5  Ndufb3  Arl3    Acsm3   Ugt2b17
Amt Ndufa10 Cstf3   Shmt2   Acy3    Gstp1
Ndufs1  Bcr Mtco2   Tpt1    Adk Pdk1
Acan    Mpc1    Bcap31  Me1 Ahcyl1  Bgn
Wasl    Cfdp1   Ndufaf2 Abcg2   Aif1l   Col4a2
Rbm3    Rps15   Gfap    H1-1    Akr1b1  Ttn
Mtnd5   Snu13   Atp6v0c Pcbp2   Amacr   Cox6b1
Ddc Lin7c   Tmed10  Ddx17   Atp2b1  Mup2
Tmc7    Ndufb5  S100a5  Slc25a1 Ca5b    Cyp2d9
Erc1    Ampd3   Ap3b1   H1-3    Cbr1    Cyp2b10
Nup35   Ipo5    Gja1    Cstb    Cbr3    Ftl1
Ndufb2  Ndufb10 Ndufb7  Tst Ces1f   Cyp8b1
Ndufb1  Ebp Ciapin1 Vat1    Cnpy2   Rbm3
Rexo2   Rabac1  Ndufb4  Glul    Col1a1  Ehd1
Cbln3   Wdfy2   Adk Rpl29   Col1a2  Cyp4a10
Arfgef3 Clcn6   Gfpt1   Pc  Col4a1  Ndufaf2
Dagla   Vbp1    Rpl28   Acy3    Col4a2  Cops2
Slc25a31    Nae1    Dgkz    Akr1a1  Cox7a1  Cyp2a4
Tph2    Gucy1a1 Tpbg    Pecr    Cpne5   Rbp1
Actbl2  Ap1s1   Gprc5b  Hnrnpd  Creg1   S100a13
Ndufb9  Ndufs1  Acyp2   Serpina1c   Crot    Gsta4
Timm10b Ndufa3  Kirrel2 Hsd17b11    Cryab   Ldhb
Syt12   Ndufs4  Ndufb11 Rpn1    Cyp2d26 Dhtkd1
Unc13c  Dnajb4  Ndufb10 Acly    Cyp2d9  Ndufb7
Serpina1d   Chmp4b  Wdr6    Atp6v1b2    Cyp2j6  Acot13
Yme1l1  Pdcd6   Mpst    Pxmp2   Cyp4a10 Fmo3
Aprt    Ndufaf2 Smdt1   Arpc3   Cyp4a14 Clta
Mtnd1   Igsf21  Sec61g  Kpna6   Cyp4b1  Eftud2
Jpt2    Stxbp3  C1qb    Atp6v1g1    Cyp7b1  Tex15
Ndufv3  Ndufb1  Ndufa2  Slc25a31    Ddc Ndufs1
Chordc1 Tpd52l2 Hddc2   Aldh3a2 Dpyd    Sult1c2
Cbln1   Sh3gl1  Agfg1   Gng12   Echdc3  Srl
Rbbp4   Pnn Manf    Tmpo    Eea1    Aimp1
Ndufaf2 Prodh   Kcnj10  Acox1   Ehhadh  Cyp2b9
Egfr    Exoc7   Ptgds   Rhoa    Folh1   Cyp4a14
Ca4 Ubxn1   Rpa3    Mapre1  Fth1    Csad
Vps33b  Nfu1    Eif3e   Gk  Ftl1    Gck
Ndufb11 Etnppl  Isg15   Rps2    Gsta2   Ss18l2
Gatad2b Spock1  Tmsb10  Atp6v1e1    Gsta3   Hebp1
Scn8a   Tagln2  Chmp1a  Sult1d1 Gsta4   Mup20
Ensa    Cbln1   Mobp    Plin1   Hacl1   Nudt7
Hkdc1   Sart3   Mtmr9   Idh1    Hmgb2   Tmem19
Ndufa5  Ndufs6  Ptprf   Cpne1   Hmgcs2  Mup18
Prpf6   Tagln   Hmox2   Prodh   Hmgn2   Selenbp1
Glipr2  Ndufs8  Ptprd   Asl Hpd Mpp1
Serpina3k   Drc7    Pmm1    Stmn1   Hsd3b2  Gstt2
Serpina1a   Hspg2   Pde6d   Abhd14b Iqgap2  Tmco1
Camk2g  Aamdc   Ipo7    Cryz    Krt18   Slc30a5
Nmnat1  Rpl26   Ndufc2  Amacr   Krt7    Hsd17b6
Tbcd    Nmral1  Csde1   Hmgb2   Krt8    Dusp3
Ckm Uck1    Itpa    Sfxn1   Ldhd    Aass
Mtif2   Lysmd2  Ndufa5  Acad11  Lmnb2   Sh3bgrl3
Hmgn5   Sars2   Actr3b  Alox12  Lsm6    Tex44
Tmem245 Cep170b Meis2   H2ac20  Mapre1  Cyp7b1
Cbx3    Bzw1    Bcas3   Csad    Mbp Tymp
Gjc3    Apex2   Ensa    Ucp1    Mpv17l  Fth1
Ugt8    Nucks1  Serpina1a   Gpd1    Mrpl22  Serpina3m
Lsmem1  Ndufb6  Serpina3k   Xab2    Mrpl23  Ces3a
Rbsn    Cp  Tnks1bp1    Adh1    Mrpl40  Trmt112
Cpped1  Npm1    Prnp    Dcps    Mrps26  Rdh16
Arpp19  Rhoa    Uap1l1  Dhx9    Myh6    Ndufs4
Ndufa6  Ndufa8  Ndufs3  Ezr Ndufa10 Hspa1b
Ndufs2  Ndufb2  Ndufv2  Sri Ndufa12 Ndufv1
Mug1    Ndufv2  Ndufs4  Ptrh2   Ndufa2  Sult1d1
Ndufs8  Nek9    Cltb    Mecr    Ndufa5  Acot3
Ndufb7  Ndufa12 Ndufa12 Sqor    Ndufa6  Ndufa12
Ndufv1  Ndufs5  Wdfy2   Echdc2  Ndufa7  Upb1
Ndufs3  Ndufa2  Ndufa10 Fabp1   Ndufa9  Ttc39c
Dnah8   Lin7b   Mtnd1   Atp2a3  Ndufaf2 Sult1a1
Ndufc2  Nek7    Ctsa    Gsr Ndufb7  Sqor
Ndufa10 Ndufs3  Ndufv3  Adrm1   Ndufs1  Aox3
Ndufs7  Ndufv1  Ndufs6  Egfr    Ndufs2  Sult2a8
Ndufs5  Ndufb4  Clybl   Gmpr    Ndufs3  Hyou1
Gcfc2   Ndufa9  Stat1   Gdap2   Ndufs4  Egfr
Ndufb3  Ndufa13 Rpl23   Ywhaq   Ndufs6  Ndufa7
Ndufb4  Ndufv3  Ndufv1  Serpina1e   Ndufv1  Ndufa6
Ndufa12 Ndufs7  Rif1    Lta4h   Ndufv2  Ndufa9
Ndufb8  Ndufb3  Ndufb8  Kifbp   Ndufv3  Srd5a1
Serpina1e   Pgam2   Capg    Psma4   Nqo1    Mpc1
Ndufb10 Ndufc2  Amt Mtr Nucks1  Keg1
Ndufs6  Ndufa6  Ndufs1  Cps1    Nudt19  Serpina3k
Fam136a Cltb    Cep170b Ehhadh  Pals2   Gsta2
Ndufv2  Ndufb11 Ndufa13 Flii    Pecr    Gstt3
Ndufa11 Pcdhga4 Vwa1    Aldob   Ppid    Metap1d
Ndufa1  Mif Tinagl1 H1-4    Ptprk   Igf2bp3
Ndufb6  Ndufb9  Ndufa11 Kyat3   Pyroxd2 Ndufs5
S100b   H1-0    Actl6b  Coasy   Rbm3    Nudt19
Acta1   Rpl29   Ndufa8  Snrpe   Reep6   Atp6v1g1
Ndufa2  Ndufs2  Lsmem1  Vta1    Retsat  Maoa
Nfia    Acad8   Ilk Osbpl8  Sephs2  Mgst3
Ndufa7  Ndufa11 Ndufa7  Agpat2  Serpina3m   Comt
Ndufa8  Ech1    Ndufs5  Clic1   Sh3bgrl3    Mat1a
Ndufs4  Chmp2a  Msn Macroh2a1   Slc12a3 Acot4
Alas1   Micos10 Nup85   Sars1   Slc22a27    Fuom
Phyhd1  S100b   Nucks1  Pabpc1  Slc22a6 Mrpl14
Mtnd4   Synpr   Tceal3  H2ac21  Slc6a19 Serpinb1a
Ndufa9  Ndufa7  Ndufa9  Nudt8   Stat1   Ces3b
Cdh10   Gnpda2  Cep152  Krt18   Tor1aip1    Prelp
Ring1   Rbfox2  Etnppl  Fuom    Ugt1a2  Mup3
Tmsb10  S1pr1   Acta1   Cyb5b       Cyp2c37
Mrpl43  mt-Nd2  Myh11   Slc25a35        Col12a1
Icmt    Rnps1   Ddx46   Rbm3        Endov
Ncln    Ap3s1   Slc44a2 Gpd1l       S100a9
Gh1 Rbx1    Ndufa6  Lrrc59      Slc25a4
Bola3   Cat Ifit1   Hmgn2       Lum
Ahsg    Atp5if1 Vps52   Tkt     Myl1
Pdlim1  Rpl18a  Rac2    Hmgb1       Col1a1
Ca3 Cyp51a1 Tmco5b  H1-2        Tpm2
Crim1   Pkn1    Bclaf1  H4c1        Myh7
Supt4h1a    Ppp3r1  Ube2v2  Stat1       Dcn
S100a9  Sltm    Atp5if1 Aldh9a1     Cyfip1
Shisa6  Chmp7   Mcm6    Elavl1      Phf11
Adcy5   Bcl7c   Vac14   H2bc7       Cse1l
Slc4a2  Lama5   Aqp1    H1-5        Col1a2
Rgs8    Prnp    Hnrnpul1    Fah     Fabp3
Asph    Nup85   Dnajc5  Aldh1l2     Acta1
Hsd17b2 Iah1    Fam136a H3-5        Myh8
Prx Mbnl2   Ublcp1  H2ac4       Ogn
Bzw2    Rbm3    Cmtm5   Ldhd        Soat2
Mboat4  Hdgfl3  Nipa1   Cryl1       Sf3b5
Cep76   Hacd2   Bzw1    Ncl     Btk
Cacna1e Asph    Tgm2    Gpd2        Cela1
Ldhd    Lpcat4  Akr1e2  Ndrg1       Uap1
    Iqsec3  Nek7    Acacb       Actn2
    Kyat3   Gabbr1  Rcc2        Tnnt3
    Hnrnpul1    Pld2    Fgb     Upp2
    Atp5f1e C1qa    Ces1f       Tgtp2
    Mrps21  Irgm1   Gclc        Tpm1
    Itpa    Try5    Pvalb       Myh1
    Sec61g  Alg2    Bdh1        Tex264
    Acads   Pik3c3  Sord        Amy2a5
    Poldip2 Dynll2  Ssb     Cela2a
    Rpl36a  Ckm Fabp4       Ckmt2
    Slc6a5  Gstm5   Apoc4       Fmod
    Fn1 Urod    Bpnt1       Des
    Usp19   Rbbp7   Fasn        Ckm
    Prtn3   Ociad1  Hcfc1       Myl3
    Sh3bgrl Rpl29   Cyp4b1      Try5
    Gclm    Ngef    Ganab       Nrde2
    Egfr    Cryab   Tkfc        Ctrb1
    Mobp    Rnps1   Acaca       Ldb3
    Abhd3   Ugdh    Cbr2
    Cirbp   Iglon5  Lipe
    Mast1   C1qc    Hnrnpa3
    Slc44a2 Eif2b4  Septin9
    Magi2   Ythdf3  Hnrnpf
    Spag7   Rps15   Hnrnpab
    C1qbp   Atg4b   Sfxn3
    Pllp    Pvalb   Ctsz
    Smad3   Lig1    Coro1c
    Exoc6b  Ehd1    S100a11
    Chtop   Fn3krp  Ssrp1
    Hmgn2   Cxadr   Ptbp1
    Abhd6   Tmem186 Rac2
    Cers1   Lsm4    Nup85
    Dsg2    Rpl10   Nherf1
    Arl8b   Ca1 Coro1a
    Rbck1   Mindy2  Atp6v1a
    Abhd4   Rcc1    Arpc1b
    Frrs1l  Sart1   Krt8
    Rnf214  H3-5    Ruvbl1
    Cdc42ep1    H1-0    Pck1
    Sf3a3   Clca3a1 Hnrnpc
    Metap2  Lpcat4  Rida
    Eif2b4  H2-K1   Ces1d
    Inpp5f  S100a13 Smarcc2
    Urod    Trappc4 Ewsr1
    Nmt2    Ca4 Akr1c21
    Alg11   S100a4  Serpina3m
    Rad23b  Hbb-b2  Chmp1a
    Drap1   Celf1   Actbl2
    Hddc2   Rap2a   Ndufa1
    App Ca8 Ndufs4
    Cpsf6   Fgg Rala
    Ppp2r5c Ube2fb  C8a
    Kiss1   Eif3k   Mug1
    Rgs6    Ctnnbl1 Rrp12
    Txndc9  Timm10  mt-Nd2
    Ehd4    Uvrag   Apoa2
    Slc25a31    Arid1a  Ndufaf2
    Cpped1  Krt8    Mtx1
    Rbm27   Bpifb9a N6amt1
    Cmtm5   Col1a2  Mrpl49
    Gh1 Actbl2  Ndufs6
    Elp3    Evl Cisd1
    Atp6v0c Rpl39   Ndufa12
    Ints3   Col6a1  Serpina3k
    Tmem68  Map4k1
    Tiam1   Atp6v1e2
    Cox7a2l Col1a1
    Arfgef3 Rpl35
    Higd1a  Cpsf6
    Ptgds   Lum
    Eny2
    Tmsb10
    Per2
    Tatdn1
    Rpl30
    Irag2
    Otud6b
    Myl1
    Akap10
    Cav1
    Fbxl20
    Adh1
    Ehhadh
    Atp4a
    Naga
    Tceal3
    Col1a1
    Celf1
    Cops6
    Vps11
    Dync1i2
    Hm13
    Eif2s2
    Fgg
    Ralgapa1
    Ipo7
    Tmod3
    Fgb
    Phpt1
    Eif3j2
    Lars2
    Snrpb2
    Mug1
    Slc4a1
    Arg1
    Celf2
    Actbl2
    Krt6a
    Fga
    Ufd1
    Fam136a
    Ntmt1
    Opalin
    Omp
    Acsf3
    Lig3
    Xpo5
    Tmem240
    Aldob
    Adprs
    Stx1a
    H2bc14
    Serpina1d
    Cldn11
    Cyp2j5
    Slc6a11
    Cpm
    Acp2
    Ca1
    Cd44
    Cplx2")
fld <- fromList(as.list(d))
upset(fld)

As shown in the figure below, only five tissues could be seen, with the exception of the kidney, even in the presence of overlapping genes between the kidney and other tissues.
.
Can you kindly assist in this regard, I do not know if there is something wrong with my input data or the script itself.

No idea how your code looks like before pasting into the web form, but in general it is not a great idea pasting +260 rows/lines of data in your script.
From what I can see in the form you have columns of data separated by different number of spaces even in the header (2 and 4 spaces).
Because of this it is impossible to figure out to which column is which.
The easiest is to keep the data in the TSV or CSV format file and make sure there are no issues with column delimiters.
EDIT
I have saved the Google Spreadsheet as de_genes.csv, edited just in case Olfactory Bulb to Olfactory_Bulb then did:
in_fn <- "de_genes.csv"

df = fread(in_fn)

de_lists <- fromList(as.list(df))

upset(de_lists, sets=c( "Brainstem",    "Cerebellum",   "Olfactory_Bulb",   "Heart", "Kidney", "Liver") , order.by = "freq")

I got:

I tried using that strategy (i.e., keeping data in csv instead of pasting it) but I am struggling with reading my csv file using the "data.table" function.        See the original data here (https://docs.google.com/spreadsheets/d/17KuGuzBjBhaRBoDCq-3SvnklBZQF1C-HcTpss1V0zKA/edit?usp=sharing)

Updated the answer. It should work now

It is indeed working, thank you.

@S'boneloGlenKhumato  could you upvote this answer and mark it as accepted please?


text:
Aligning FASTQs to FASTA reference
I'd like to align some FASTQs to an average mtDNA FASTA file that I have downloaded so I can have the human mtDNA isolated from those FASTQs. For that, I used bowtie2. Can I expect that after running the command bowtie2 -x ref_index -U seq.fastq.gz -S seq.sam, my seq.sam will only contain the DNA that I'm interested in (i.e., Human mtDNA)?
I have this example below from a class in which is supposed that the human DNA is discarded (Opposite of what I want). But I don't understand the roles of the line 3. Is it necessary to filter out undesired DNA?

bowtie2 -p $THREADS -x host_DB -1 $FW_READS -2 $RV_READS --un-conc-gz SAMPLE_host_removed > $STUDY-data.sam

samtools view -bS $STUDY-data.sam > $STUDY-data.bam

samtools view -b -f 12 -F 256 $STUDY-data.bam > $STUDY-data-no-human.bam

samtools sort -n -o $STUDY-sorted-data.bam $STUDY-data-no-human.bam 1>> $REPORTS_DIR/report_stdout.txt

samtools fastq -1 $STUDY-reads-r1.fq -2 $STUDY-reads-r2.fq $STUDY-sorted-data.bam 1>> $REPORTS_DIR/report_stdout.txt

Can I expect that after running the command bowtie2 -x ref_index -U seq.fastq.gz -S seq.sam, my seq.sam will only contain the DNA that I'm interested in (i.e., Human mtDNA)?

Not quite as that depends on what you had in your input. If you had a lot of non-human mitochondrial sequences, for example, those might still be there if they align well enough against your reference. Essentially, what you can assume here is that the SAM file will have those input sequences that aligned well against your mt DNA reference. Whether that will be a clean set of human mitochondrial sequences depends entirely on what you had as input going in.
However, as a general rule, yes, aligning against a specific target should leave you with the reads that align to that target which in your case means human mitochondrial sequences.

And here is a brief explanation of what each command does:

bowtie2 -p $THREADS -x host_DB -1 $FW_READS -2 $RV_READS --un-conc-gz SAMPLE_host_removed > $STUDY-data.sam
This will align the input fastq sequences against the reference to produce a SAM file.

samtools view -bS $STUDY-data.sam > $STUDY-data.bam
This simply converts the SAM file to a BAM file. Essentially a compressed SAM.

samtools view -b -f 12 -F 256 $STUDY-data.bam > $STUDY-data-no-human.bam
Here, you are filtering your bam. The -b flag just means the output will be BAM, and the -f and -F control what alignments are kept or discarded respectively depending on the SAM flag. This... is a bit complicated, you need to read the SAM specification to understand fully, but briefly, each alignment in a SAM/BAM file has a binary flag whose value gives you information about the alignment. You can use the helpful "Explain SAM flags" page from the Broad Institute to understand what each number means.
12 means "read unmapped and mate unmapped" so -f 12 means "keep only alignments where the bits given are set". this will match read pairs where neither pair has been mapped. Next, -F is the inverse option and means "remove alignments where the bits given are set". 256 means "this is not a primary alignment". Taken together then, this command would only keep unmapped reads, excluding reads that are unmapped as a secondary alignment. In other words, it will result in a bam file with only those reads that couldn't be mapped at all. Presumably, in the example you found this from, this is to collect the non-human reads, i.e. those that don't match the reference genome.

samtools sort -n -o $STUDY-sorted-data.bam $STUDY-data-no-human.bam 1>> $REPORTS_DIR/report_stdout.txt
This just sorts the bam file. It doesn't remove anything.

samtools fastq -1 $STUDY-reads-r1.fq -2 $STUDY-reads-r2.fq $STUDY-sorted-data.bam 1>> $REPORTS_DIR/report_stdout.txt
And this just extracts fastq files from the BAM alignment.


text:
Tools and software to visualize the moving cell compartments in 3D
Also asked on biostars
I am looking for tools and software that could help me visualize moving cell compartments in 3D where I can go from the cell membrane to the ER, then out with the protein into the matrix. In short a 3D visualization of the cell environment.
Since I am a newbie, I would appreciate any tool/software suggestion for me to try.

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.


text:
How to extract mutations in a particular gene from the VCF file of whole genome
Also posted on biostars
I would like to extract all the mutations in a single gene of interest from the VCF file. For example I would like to extract all the mutations in the NFKB1 gene from the whole genome VCF file.

What have you tried? bcftools view has an option to extract by coordinates or by annotation if your VCF had the gene annotation. If you used custom annotation tools, they may have their own filtering utilities.

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.

Find the location of the gene of interest and then filter the vcf for that location using vcftools, bcftools, plink, or some other software. There is a lot of good documentation available online for these tools.

As mentioned, there are many options to go (eg. vcftools and bcftools). The Broad Institute provides a valuable toolset (GATK) for this very purpose and extended procedures for various variant selection tasks.
SelectVariants function can be good for this purpose, though the location of a gene of interest should be addressed:
https://gatk.broadinstitute.org/hc/en-us/articles/13832694334235-SelectVariants


text:
WGCNA: Extract genes within the different modules to create a dataframe with two columns
My inquiry pertains to WGCNA. After identifying the modules, I would like to extract genes within the different modules to create a dataframe with two columns: Gene_a and Gene_b. The values in these columns should represent co-expressed genes. An illustrative example is as follows:

In this example, the gene XX62730x is co-expressed with XCGR1Y. And YYYYY is co-expressed with ZZZZZZ.
However, I am facing difficulties as the official WGCNA page is no longer accessible.
I would greatly appreciate any help you can provide. Thank you in advance!
Best,
Nana

The archive.org copy might help while you're trying to figure out the answer: https://web.archive.org/web/20231002045159/https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/

@RamRS thanks for the link!

Unfortunately, the links on that link  do not work.

I tested a couple and they worked - you need to learn to navigate archive.org to be able to leverage those links the best way.

U can follow WGCNA R packages

You can get this information, more or less, using the functions exportNetworkToVisANT and exportNetworkToCytoscape. You will find some extra columns that you can ignore or delete from the resulting table. We are still working with UCLA to restore the WGCNA page including the tutorials, but the tutorials can temporarily be accessed using this dropbox link.
To use either one of the two functions, you will need the network matrix for which we usually use TOM. You will probably want to set a threshold for connections to be included, otherwise you will get a fully connected network and the resulting table will (usually) be huge.
Hope this helps!

thank you very much, that really helped me


text:
Bedtools genomecov gives me coverage files showing zeros at all nucleotide positions
Also posted on biostars
I have bam files that I've obtained through STAR alignment. Here is the first 10 lines of one of my bam files:
GWNJ-1013:671:GW2304233104th:2:1126:15501:17879 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1126:15130:18427 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1206:9272:29857  0       chromosome      1       255     28S108M14S      *       0       0       GGCGGTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGAGATCGGAAGAGCG  FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:106        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:30734:10692 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:31412:33411 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFF::FFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFF:FFFF,FFFFFFF:FFFFFFFFF:FFFFFFFFFFFFFFF:FFF::FFFFFFFF,FFFFFFFF:FFFFFFFFFFF:F,FFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1330:21938:10708 0       chromosome      1       255     16S134M *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF,FFFFFFFFFFFFFFFFF:F  NH:i:1  HI:i:1  AS:i:132        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:17074:4726  0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:11939:11584 0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1372:24740:30843 0       chromosome      1       255     16S48M  *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAA        FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        NH:i:1  HI:i:1  AS:i:47 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1401:11107:19319 0       chromosome      1       255     5S144M1S        *       0       0       ATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGGAATTGA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF  NH:i:1  HI:i:1  AS:i:142        nM:i:0

This is the script I'm using that contains the genomecov -d -pc -ibam command:
BAM_DIR='./sorted_bam_files'

# Loop over each BAM file in the directory
for file in "$BAM_DIR"/*_high_mapq_reads.sorted.bam
do
    echo "Calculating coverage for $file..."
    # Get the basename of the file
    base=$(basename "$file")
    # Calculate coverage
    bedtools genomecov -d -pc -ibam $file > "${base%.sorted.bam}.coverage"
done

This is the output I keep getting (first 10 lines only shown). I've tried this with at least 30 other bam files but the output is always the same.
chromosome      1       0
chromosome      2       0
chromosome      3       0
chromosome      4       0
chromosome      5       0
chromosome      6       0
chromosome      7       0
chromosome      8       0
chromosome      9       0
chromosome      10      0

Any ideas on why this keeps happening? I don't think there's something wrong with my bam files because I was able to use them for DGE analyses and also checked to see if most reads were getting unmapped but they weren't. Any help is appreciated.

Why is your first chromosome called "chromosome", that's odd.

@ATpoint also second column FLAG == 0 and the fifth MAPQ ==255 meaning: "A value 255 indicates that the mapping quality is not available.". Then the CIGARs looking like the read was indeed mapped 11S82M etc.

@ATpoint The data I have is from bacteria, so I'm only working with 1 chromosome

option -pc is for
-pc     Calculate coverage of pair-end fragments.

but your reads are single-end reads

What tells you they're single ends? Is it because you know this would happen with -pc and single end files or was it something in the bam?

"What tells you they're single ends? Is it b" : the sam flags of each alignemt is zero. Paired would have the '1' bit set . see https://broadinstitute.github.io/picard/explain-flags.html

Using that page, I can see that 11S139M, the first one, is "read paired", "read mapped in proper pair" and "mate unmapped". Could you explain a bit more please? Ideally, adding the detail in the answer so the next user who is confused can also find the info?

the second column is the read flag (section 1.4 of https://samtools.github.io/hts-specs/SAMv1.pdf) . it's always 0  in your example.  0 is single end read mapped on plus strand.

Ah thanks, yes, I was looking at the cigar string, duh. It would be great if you could add some more detail to the answer to make all this clear. Especially since I don't see mention of 0 meaning "single end" either in the specs (but perhaps I missed it) or when using the explain-flags page (which just ignores 0).

"I don't see mention of 0 meaning "single end" either in the specs (but perhaps I missed it)" it's a array of bits https://en.wikipedia.org/wiki/Bit_array . The first bit is for "Read Paired". As it's not set, the read is not part of a "Read Paired", hence it's single end.

@Pierre Some of my bam files have a mix of 0 and 16 for my read flags. So to my understanding I have a mix of single-end and paired-end reverse strand reads. Do I include or omit -pc?


text:
How to shade text in a table
I have currently made a tool that takes input sequences and maps them to a reference sequence highlighting it. You can then save the word document and it creates a table of the sequence name and the sequences.
What I want is to

Highlight the sequences in specific colours.
Highlight those sequences in the reference sequence.

Below is the create a word function.
def create_word_document(output_text, reference_sequence):
    document = Document()
    reverse_complemented_sequences = {}

    # Add reference sequence
    document.add_heading("Reference Sequence", level=1)
    document.add_paragraph(reference_sequence, style='BodyText')

    # Add a table for sequences
    table = document.add_table(rows=1, cols=2, style='Table Grid')
    table.autofit = True
    table.cell(0, 0).text = "Sequence Name"
    table.cell(0, 1).text = "Highlighted Sequence"

    # Colors for highlighting sequences
    highlight_colors = {
        "F3": "FF0000",  # Red for F3
        "B3": "0000FF",  # Blue for B3
        "F2": "FFC0CB", # Pink for F2
        "F1c": "0000FF", # Blue for F1c
        "B2": "FFFF00", # Yellow for B2
        "B1c": "00008B", # Dark blue for B1c
        "LF": "FFA500", # Orange for LF
        "LB": "800080", # Purple for LB
        # Change colors as needed
    }

    for seq_name, highlighted_sequence in output_text.highlighted_sequences.items():
        # Use the reverse complemented sequence for this section
        reversed_sequence = reverse_complement_if_needed(seq_name, highlighted_sequence)
        reverse_complemented_sequences[seq_name] = reversed_sequence

        # Add sequence information to the table
        row_cells = table.add_row().cells
        row_cells[0].text = f"{seq_name} Sequence"
        row_cells[1].text = reversed_sequence

        # Add shading to the cell containing the sequence text if available
        if seq_name in highlight_colors:
            color_code = highlight_colors[seq_name]
            apply_shading_to_cell(row_cells[1], reversed_sequence, color_code)

        # Add FIP and BIP sequences after B3
        if seq_name == "B3":
            f1c_sequence = reverse_complement_if_needed("F1c", output_text.highlighted_sequences["F1c"])
            f2_sequence = output_text.highlighted_sequences["F2"]
            bip_sequence = reverse_complement_if_needed("B1c", output_text.highlighted_sequences["B1c"])
            b2_sequence = output_text.highlighted_sequences["B2"]

            # Add a row for FIP sequence
            row_cells_fip = table.add_row().cells
            row_cells_fip[0].text = "FIP Sequence"
            row_cells_fip[1].text = f"{f1c_sequence}-{f2_sequence}"
            # Apply shading to FIP sequence if available
            if "F1c" in highlight_colors:
                apply_shading_to_cell(row_cells_fip[1], row_cells_fip[1].text, highlight_colors["F1c"])

            # Add a row for BIP sequence
            row_cells_bip = table.add_row().cells
            row_cells_bip[0].text = "BIP Sequence"
            row_cells_bip[1].text = f"{bip_sequence}-{b2_sequence}"
            # Apply shading to BIP sequence if available
            if "B1c" in highlight_colors:
                apply_shading_to_cell(row_cells_bip[1], row_cells_bip[1].text, highlight_colors["B1c"])

    # Save the document
    document.save("LAMP Primer mapping.docx")

def apply_shading_to_cell(cell, text, color_code):
    # Clear existing text
    cell.text = ''
    # Add the text to the cell
    paragraph = cell.paragraphs[0]
    run = paragraph.add_run(text)
    # Create shading element
    shading_elm = OxmlElement('w:shd')
    shading_elm.set(qn('w:val'), 'clear')
    shading_elm.set(qn('w:color'), 'auto')
    shading_elm.set(qn('w:fill'), color_code)  # Set fill color
    # Append shading element to run properties
    run._r.append(shading_elm)
    # Set font size for the entire paragraph
    paragraph.style.font.size = Pt(12)

Full code:
import tkinter as tk
from Bio.Seq import Seq
from Bio import pairwise2
from docx import Document
from docx.shared import Inches
from docx.shared import Pt
from docx.oxml.shared import OxmlElement
from docx.oxml.ns import qn

class HighlightedOutputText(tk.Text):
    def __init__(self, *args, **kwargs):
        tk.Text.__init__(self, *args, **kwargs)
        self.highlighted_sequences = {}

    def map_and_highlight(self, seq_name, reference_sequence, input_sequence):
        best_alignment = self.find_best_alignment(reference_sequence, input_sequence)

        if best_alignment is not None:
            reference_aligned, input_aligned, start_position, end_position = best_alignment

            # Configure highlight tag based on sequence name
            highlight_tag = f"{seq_name.lower()}_highlight"
            self.tag_configure(highlight_tag, background=get_highlight_color(seq_name))

            # Apply the color tag to the aligned region in the reference sequence
            self.insert(tk.END, f"\n{seq_name} Input sequence:\n{input_sequence}\n\n")
            self.insert(tk.END, "\nHighlighted sequence:\n")
            self.insert(tk.END, reference_sequence[:start_position], "default")
            self.insert(tk.END, reference_sequence[start_position:end_position], highlight_tag)
            self.insert(tk.END, reference_sequence[end_position:], "default")

            # Save highlighted sequence for later use
            self.highlighted_sequences[seq_name] = reference_sequence[start_position:end_position]
        else:
            self.insert(tk.END, f"\n{seq_name} Input sequence:\n{input_sequence}\n\n")
            self.insert(tk.END, "\nNo matching region found in the reference sequence.\n", "default")

    def find_best_alignment(self, reference_sequence, input_sequence):
        alignments_forward = pairwise2.align.localxx(reference_sequence, input_sequence)
        alignments_reverse = pairwise2.align.localxx(reference_sequence, str(Seq(input_sequence).reverse_complement()))

        if not alignments_forward and not alignments_reverse:
            print(f"No alignment found for input sequence: {input_sequence}")
            return None

        # Choose the best alignment based on the score
        best_alignment = max(alignments_forward + alignments_reverse, key=lambda a: a.score)

        # Extract aligned sequences
        reference_aligned = best_alignment[0]
        input_aligned = best_alignment[1]

        # Find the start and end positions of the input sequence in the reference
        start_position = reference_aligned.find(input_aligned.replace('-', '').lstrip())
        end_position = start_position + len(input_aligned.replace('-', '').rstrip())

        return reference_aligned, input_aligned, start_position, end_position

    def clear_output(self):
        self.delete(1.0, tk.END)
        self.highlighted_sequences = {}

def reverse_complement_if_needed(seq_name, input_sequence):
    # Check if the sequence name indicates reverse complementing is needed
    if seq_name in ["B3", "F1c", "B2", "LF"]:
        return str(Seq(input_sequence).reverse_complement())
    else:
        return input_sequence

def get_highlight_color(seq_name):
    # Return RGB values based on the sequence name
    if seq_name == "F3":
        return "#808080"  # Gray
    elif seq_name == "B3":
        return "#404040"  # Dark Gray
    elif seq_name == "LF":
        return "#008000"  # Green
    elif seq_name == "LB":
        return "#FF0000"  # Red
    elif seq_name == "F2":
        return "#40E0D0"  # Turquoise
    elif seq_name == "F1c":
        return "#FFFF00"  # Yellow
    elif seq_name == "B2":
        return "#FFC0CB"  # Pink
    elif seq_name == "B1c":
        return "#00008B"  # Dark Blue
    else:
        return "#FFA500"  # Default color for other sequences

def create_word_document(output_text, reference_sequence):
    document = Document()
    reverse_complemented_sequences = {}

    # Add reference sequence
    document.add_heading("Reference Sequence", level=1)
    document.add_paragraph(reference_sequence, style='BodyText')

    # Add a table for sequences
    table = document.add_table(rows=1, cols=2, style='Table Grid')
    table.autofit = True
    table.cell(0, 0).text = "Sequence Name"
    table.cell(0, 1).text = "Highlighted Sequence"

    # Colors for highlighting sequences
    highlight_colors = {
        "F3": "FF0000",  # Red for F3
        "B3": "0000FF",  # Blue for B3
        "F2": "FFC0CB", # Pink for F2
        "F1c": "0000FF", # Blue for F1c
        "B2": "FFFF00", # Yellow for B2
        "B1c": "00008B", # Dark blue for B1c
        "LF": "FFA500", # Orange for LF
        "LB": "800080", # Purple for LB
        # Change colors as needed
    }

    for seq_name, highlighted_sequence in output_text.highlighted_sequences.items():
        # Use the reverse complemented sequence for this section
        reversed_sequence = reverse_complement_if_needed(seq_name, highlighted_sequence)
        reverse_complemented_sequences[seq_name] = reversed_sequence

        # Add sequence information to the table
        row_cells = table.add_row().cells
        row_cells[0].text = f"{seq_name} Sequence"
        row_cells[1].text = reversed_sequence

        # Add shading to the cell containing the sequence text if available
        if seq_name in highlight_colors:
            color_code = highlight_colors[seq_name]
            apply_shading_to_cell(row_cells[1], reversed_sequence, color_code)

        # Add FIP and BIP sequences after B3
        if seq_name == "B3":
            f1c_sequence = reverse_complement_if_needed("F1c", output_text.highlighted_sequences["F1c"])
            f2_sequence = output_text.highlighted_sequences["F2"]
            bip_sequence = reverse_complement_if_needed("B1c", output_text.highlighted_sequences["B1c"])
            b2_sequence = output_text.highlighted_sequences["B2"]

            # Add a row for FIP sequence
            row_cells_fip = table.add_row().cells
            row_cells_fip[0].text = "FIP Sequence"
            row_cells_fip[1].text = f"{f1c_sequence}-{f2_sequence}"
            # Apply shading to FIP sequence if available
            if "F1c" in highlight_colors:
                apply_shading_to_cell(row_cells_fip[1], row_cells_fip[1].text, highlight_colors["F1c"])

            # Add a row for BIP sequence
            row_cells_bip = table.add_row().cells
            row_cells_bip[0].text = "BIP Sequence"
            row_cells_bip[1].text = f"{bip_sequence}-{b2_sequence}"
            # Apply shading to BIP sequence if available
            if "B1c" in highlight_colors:
                apply_shading_to_cell(row_cells_bip[1], row_cells_bip[1].text, highlight_colors["B1c"])

    # Save the document
    document.save("LAMP Primer mapping.docx")

def apply_shading_to_cell(cell, text, color_code):
    # Clear existing text
    cell.text = ''
    # Add the text to the cell
    paragraph = cell.paragraphs[0]
    run = paragraph.add_run(text)
    # Create shading element
    shading_elm = OxmlElement('w:shd')
    shading_elm.set(qn('w:val'), 'clear')
    shading_elm.set(qn('w:color'), 'auto')
    shading_elm.set(qn('w:fill'), color_code)  # Set fill color
    # Append shading element to run properties
    run._r.append(shading_elm)
    # Set font size for the entire paragraph
    paragraph.style.font.size = Pt(12)

# GUI setup
root = tk.Tk()
root.title("Sequence Alignment")

# Entry widgets for reference sequence
label_reference = tk.Label(root, text="Reference Sequence:")
entry_reference = tk.Entry(root, width=30)

# Entry widgets for multiple sequences
entries = []
output_text = HighlightedOutputText(root, height=30, width=40)
output_text.configure(font=("TkDefaultFont", 10))

for seq_name in ["F3", "B3", "F2", "F1c", "B2", "B1c", "LF", "LB"]:
    label = tk.Label(root, text=f"{seq_name} Sequence:")
    entry = tk.Entry(root, width=30)
    entries.append((seq_name, entry))
    label.grid(row=len(entries), column=0, padx=10, pady=5)
    entry.grid(row=len(entries), column=1, padx=10, pady=5)

def on_align_button_click():
    reference_sequence = entry_reference.get()

    output_text.clear_output()  # Clear previous output

    for seq_name, entry in entries:
        input_sequence = entry.get()
        input_sequence = reverse_complement_if_needed(seq_name, input_sequence)
        output_text.map_and_highlight(seq_name, reference_sequence, input_sequence)

    # Display all highlighted sequences in the output window
    for seq_name, highlighted_sequence in output_text.highlighted_sequences.items():
        output_text.insert(tk.END, f"\n{seq_name} Input sequence:\n{highlighted_sequence}\n")

# Button to trigger alignment
align_button = tk.Button(root, text="Align", command=on_align_button_click)

# Button to create a Word document
document_button = tk.Button(root, text="Create Word Document", command=lambda: create_word_document(output_text,entry_reference.get()))

# Grid layout
label_reference.grid(row=0, column=0, padx=10, pady=5)
entry_reference.grid(row=0, column=1, padx=10, pady=5)

align_button.grid(row=len(entries) + 1, column=0, columnspan=2, pady=10)
document_button.grid(row=len(entries) + 2, column=0, columnspan=2, pady=10)
output_text.grid(row=len(entries) + 3, column=0, columnspan=2, padx=10, pady=5)

root.mainloop()
```

this is a pretty good question in terms of detail and complete code etc. However i think it is more of a stackoverflow question about writing word docs in python.

I agree but I have done this (years ago) - there's validity to it being here. However, reading that amount of code for customised colouring for something I can hardly remember now isn't my idea of fun.

I dunno about Word, but you can definitely do this kind of thing in Excel via Python with openpyxl, if that helps.


text:
Converting EntrezGene IDs to Ensemble IDs
I am trying to convert my gene_id accessions given by the RNAseq as EntrezGene gene numbers into Ensemble annotation. I need the Ensemble annotation for further GO analysis. I am using BiomaRT package in R studio.
In order to create a list with all my EntrezGene IDs I did the following:
gene_data <- read_excel('gene_diff_control vs treatment.xls')
mygenes <- gene_data$gene_id  # gene_id is the column containing the EntrezGene IDs

# here is an example of the gene_id`s I have in my file: gene_diff_control vs treatmen.xls"
# 27304
# 8228
# 83698
# 7283
# and so on till row number 16,620

# then I use the following code:
hsmart <- useMart(dataset = "hsapiens_gene_ensembl", biomart = "ensembl")

mapping <- getBM(
    attributes = c('ensembl_gene_id', 'ensembl_transcript_id', 'entrezgene_id', 'hgnc_symbol'),
    filters = 'hgnc_symbol',
    values = mygenes,
    mart = hsmart
)

#Eventually the mapping gets 0 entries:
[1] ensembl_gene_id       ensembl_transcript_id entrezgene_id
[4] hgnc_symbol
<0 rows> (or 0-length row.names)

# So I modified my code in order to use the 'gene_symbol' instead

mygenes <- gene_data$gene_symbol
hsmart <- useMart(dataset = "hsapiens_gene_ensembl", biomart = "ensembl")
mapping <- getBM(
    attributes = c('ensembl_gene_id', 'ensemble_gene_id_version', 'ensembl_transcript_id', 'entrezgene_id', 'hgnc_symbol'),
    filters = 'hgnc_symbol',
    values = mygenes,
    mart = hsmart
)

But then, the resulting mapping output contain more than 300,000 rows of mapping data between EntrezGene ID and Ensemble gene IDs not related to my dataset.
What am I doing wrong?

Change filter to entrez_gene_id - the IDs you have are not hgnc_symbol.

I am not 100% sure but the 300,000 rows might be corresponding to the transcripts, you have 'ensembl_transcript_id' within the attributes argument when calling getBM().

Something went wrong with this post. I had two very useful comments and now they are gone!!! If any of you guys see this post again, please write again what you think about the issue described above.

Now I`ve simplified the code to the following lines:

library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
gene_data <- read_excel('MyData.xls')
mygenes <- gene_data$gene_id

hsmart <- useMart(dataset = "hsapiens_gene_ensembl", biomart = "ensembl")
mapping <- getBM(

     attributes = 'ensembl_gene_id',

     filters = 'entrezgene_id', #

     values = mygenes,

     mart = hsmart

 )

my data set contains 16620 rows. The output of mapping is 18773 rows with corresponding ensemble gene IDs.
What is wrong?

Please update the question with the changes you have made. I have undeleted the comments. If they are indeed helpful and solve the issue, it would be great if you could post an answer explaining how you solved it with their guidance.


text:
Unable to install prokka using homebrew or conda
I have tried to install prokka using homebrew using the code brew install prokka (I also tried brew install brewsci/bio/prokka) and it gave this error:
curl: (7) Failed to connect to mbio-serv2.mbioekol.lu.se port 80 after 75009 ms: Couldn't connect to server
Error: prokka: Failed to download resource "aragorn"
Download failed: http://mbio-serv2.mbioekol.lu.se/ARAGORN/Downloads/aragorn1.2.38.tgz

I then tried installing with conda using the code
conda create -n prokka_env
conda activate prokka_env
mamba install prokka

and it gave the error:
Could not solve for environment specs
The following package could not be installed
└─ prokka does not exist (perhaps a typo or a missing channel).

Prokka's github page  github page has:
conda install -c conda-forge -c bioconda -c defaults prokka

edit
Seems that there may be another way to do it:
stepwise components install

I've tried that as well but I get the error:
'''CondaValueError: You have chosen a non-default solver backend (libmamba) but it was not recognized. Choose one of: classic'''

I have tried on Mac (
ProductName:  macOS
ProductVersion:  14.5
BuildVersion:  23F79 ).
No errors.

conda config --set restore_free_channel true

@LucyKelly  • Choose one of: classic •

libmamba vs classic -->  https://conda.github.io/conda-libmamba-solver/user-guide/libmamba-vs-classic/

@LucyKelly do you have mamba installed and on the PATH?


text:
Running multiple sequence using spades.py on an HPC
I am trying to run a batch script to assemble sequences using the spades.py command on HPC. I tested the script with one sequence (1 forward and 1 reverse read) and this worked. I have 100 sequences in total (200 files) and am trying to assemble all of them without having to resubmit the script each time. This is the original script that worked:
#!/bin/bash

#SBATCH --job-name=pangenome_spades
#SBATCH --account=lilianasalvador
#SBATCH --partition=standard
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=5gb
#SBATCH --time=01:00:00

spades.py -1 SRR8599992_1.fastq.gz -2 SRR8599992_2.fastq.gz -o /home/u32/lucykelly/spades_test -t 16 --careful

This is the script that I've been trying to run and the error it's giving me:
#!/bin/bash

#SBATCH --job-name=pangenome_spades
#SBATCH --account=lilianasalvador
#SBATCH --partition=standard
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=5gb
#SBATCH --time=02:00:00

module load python/3.9
module load spades/3.15.0-hpc

spades.py -1 *_1.fastq.gz -2 *_2.fastq.gz -o /groups/lilianasalvador/noah_100/SPAdes_output -t 16 --careful

SPAdes genome assembler v3.15.0-hpc

Usage: spades.py [options] -o <output_dir>
spades.py: error: Please specify option (e.g. -1, -2, -s, etc)) for the following paths: SRR8599995_1.fastq.gz, SRR8599997_1.fastq.gz, SRR8599998_1.fastq.gz, SRR8600006_1.fastq.gz, SRR8600007_1.fastq.gz, SRR8600008_1.fastq.gz, SRR8600014_1.fastq.gz, SRR8600015_1.fastq.gz, SRR8600020_1.fastq.gz, SRR8600022_1.fastq.gz, SRR8600026_1.fastq.gz, SRR8600031_1.fastq.gz, SRR8600044_1.fastq.gz, SRR8600045_1.fastq.gz, SRR8600048_1.fastq.gz, SRR8600049_1.fastq.gz, SRR8600052_1.fastq.gz, SRR8600055_1.fastq.gz, SRR8600058_1.fastq.gz, SRR8600059_1.fastq.gz, SRR8600064_1.fastq.gz, SRR8600070_1.fastq.gz, SRR8600072_1.fastq.gz, SRR8600075_1.fastq.gz, SRR8600079_1.fastq.gz, SRR8600092_1.fastq.gz, SRR8600093_1.fastq.gz, SRR8600098_1.fastq.gz, SRR8600100_1.fastq.gz, SRR8600101_1.fastq.gz, SRR8600102_1.fastq.gz, SRR8600103_1.fastq.gz, SRR8600105_1.fastq.gz, SRR8600110_1.fastq.gz, SRR8600113_1.fastq.gz, SRR8600117_1.fastq.gz, SRR8600121_1.fastq.gz, SRR8600122_1.fastq.gz, SRR8600123_1.fastq.gz, SRR8600128_1.fastq.gz, SRR8600133_1.fastq.gz, SRR8600134_1.fastq.gz, SRR8600136_1.fastq.gz, SRR8600145_1.fastq.gz, SRR8600147_1.fastq.gz, SRR8600160_1.fastq.gz, SRR8600164_1.fastq.gz, SRR8600166_1.fastq.gz, SRR8600167_1.fastq.gz, SRR8600174_1.fastq.gz, SRR8600182_1.fastq.gz, SRR8600184_1.fastq.gz, SRR8600185_1.fastq.gz, SRR8600188_1.fastq.gz, SRR8600189_1.fastq.gz, SRR8600190_1.fastq.gz, SRR8600193_1.fastq.gz, SRR8600194_1.fastq.gz, SRR8600216_1.fastq.gz, SRR8600222_1.fastq.gz, SRR8600223_1.fastq.gz, SRR8600226_1.fastq.gz, SRR8600228_1.fastq.gz, SRR8600235_1.fastq.gz, SRR8600236_1.fastq.gz, SRR8600243_1.fastq.gz, SRR8600244_1.fastq.gz, SRR8600251_1.fastq.gz, SRR8600252_1.fastq.gz, SRR8600258_1.fastq.gz, SRR8600259_1.fastq.gz, SRR8600264_1.fastq.gz, SRR8600267_1.fastq.gz, SRR8600268_1.fastq.gz, SRR8600270_1.fastq.gz, SRR8600272_1.fastq.gz, SRR8600274_1.fastq.gz, SRR8600277_1.fastq.gz, SRR8600278_1.fastq.gz, SRR8600281_1.fastq.gz, SRR8600282_1.fastq.gz, SRR8600283_1.fastq.gz, SRR8600285_1.fastq.gz, SRR8600287_1.fastq.gz, SRR8600288_1.fastq.gz, SRR8600291_1.fastq.gz, SRR8600294_1.fastq.gz, SRR8600296_1.fastq.gz, SRR8600302_1.fastq.gz, SRR8600308_1.fastq.gz, SRR8600309_1.fastq.gz, SRR8600311_1.fastq.gz, SRR8600319_1.fastq.gz, SRR8600327_1.fastq.gz, SRR8600328_1.fastq.gz, SRR8600329_1.fastq.gz, SRR8600331_1.fastq.gz, SRR8600339_1.fastq.gz, SRR8599995_2.fastq.gz, SRR8599997_2.fastq.gz, SRR8599998_2.fastq.gz, SRR8600006_2.fastq.gz, SRR8600007_2.fastq.gz, SRR8600008_2.fastq.gz, SRR8600014_2.fastq.gz, SRR8600015_2.fastq.gz, SRR8600020_2.fastq.gz, SRR8600022_2.fastq.gz, SRR8600026_2.fastq.gz, SRR8600031_2.fastq.gz, SRR8600044_2.fastq.gz, SRR8600045_2.fastq.gz, SRR8600048_2.fastq.gz, SRR8600049_2.fastq.gz, SRR8600052_2.fastq.gz, SRR8600055_2.fastq.gz, SRR8600058_2.fastq.gz, SRR8600059_2.fastq.gz, SRR8600064_2.fastq.gz, SRR8600070_2.fastq.gz, SRR8600072_2.fastq.gz, SRR8600075_2.fastq.gz, SRR8600079_2.fastq.gz, SRR8600092_2.fastq.gz, SRR8600093_2.fastq.gz, SRR8600098_2.fastq.gz, SRR8600100_2.fastq.gz, SRR8600101_2.fastq.gz, SRR8600102_2.fastq.gz, SRR8600103_2.fastq.gz, SRR8600105_2.fastq.gz, SRR8600110_2.fastq.gz, SRR8600113_2.fastq.gz, SRR8600117_2.fastq.gz, SRR8600121_2.fastq.gz, SRR8600122_2.fastq.gz, SRR8600123_2.fastq.gz, SRR8600128_2.fastq.gz, SRR8600133_2.fastq.gz, SRR8600134_2.fastq.gz, SRR8600136_2.fastq.gz, SRR8600145_2.fastq.gz, SRR8600147_2.fastq.gz, SRR8600160_2.fastq.gz, SRR8600164_2.fastq.gz, SRR8600166_2.fastq.gz, SRR8600167_2.fastq.gz, SRR8600174_2.fastq.gz, SRR8600182_2.fastq.gz, SRR8600184_2.fastq.gz, SRR8600185_2.fastq.gz, SRR8600188_2.fastq.gz, SRR8600189_2.fastq.gz, SRR8600190_2.fastq.gz, SRR8600193_2.fastq.gz, SRR8600194_2.fastq.gz, SRR8600216_2.fastq.gz, SRR8600222_2.fastq.gz, SRR8600223_2.fastq.gz, SRR8600226_2.fastq.gz, SRR8600228_2.fastq.gz, SRR8600235_2.fastq.gz, SRR8600236_2.fastq.gz, SRR8600243_2.fastq.gz, SRR8600244_2.fastq.gz, SRR8600251_2.fastq.gz, SRR8600252_2.fastq.gz, SRR8600258_2.fastq.gz, SRR8600259_2.fastq.gz, SRR8600264_2.fastq.gz, SRR8600267_2.fastq.gz, SRR8600268_2.fastq.gz, SRR8600270_2.fastq.gz, SRR8600272_2.fastq.gz, SRR8600274_2.fastq.gz, SRR8600277_2.fastq.gz, SRR8600278_2.fastq.gz, SRR8600281_2.fastq.gz, SRR8600282_2.fastq.gz, SRR8600283_2.fastq.gz, SRR8600285_2.fastq.gz, SRR8600287_2.fastq.gz, SRR8600288_2.fastq.gz, SRR8600291_2.fastq.gz, SRR8600294_2.fastq.gz, SRR8600296_2.fastq.gz, SRR8600302_2.fastq.gz, SRR8600308_2.fastq.gz, SRR8600309_2.fastq.gz, SRR8600311_2.fastq.gz, SRR8600319_2.fastq.gz, SRR8600327_2.fastq.gz, SRR8600328_2.fastq.gz, SRR8600329_2.fastq.gz, SRR8600331_2.fastq.gz, SRR8600339_2.fastq.gz
```

This is easy. I'll respond over the weekend.

There is brand new spades release with a bunch of improvements. You may check if it is compatible with your HPC setup. https://github.com/ablab/spades/releases/tag/v4.0.0

I have some time now, so I'm going to take a shot at this, @M__
You need to run spades per read pair. When you specify *_1.fastq.gz as the -1 param, you're giving it every _1 file at once, not one _1 file per spades command.
I'd recommend this:
for smpl in $(ls *_1.fastq.gz | xargs -I infile basename infile _1.fastq.gz)
do
    cat > "./${smpl}_spades.slurm" <<EOF
#!/bin/bash

#SBATCH --job-name=pangenome_spades.${smpl}
#SBATCH --account=lilianasalvador
#SBATCH --partition=standard
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=5gb
#SBATCH --time=02:00:00

spades.py -1 ${smpl}_1.fastq.gz -2 ${smpl}_2.fastq.gz -o groups/lilianasalvador/noah_100/SPAdes_output_${smpl} -t 16 --careful
EOF
done

Run a sanity check on a few of the .slurm files, then
ls *.slurm | xargs -I job_file sbatch job_file

Cool, thanks. I was going to solve this via a formal Slurm array job, but this will work.

Please do add an array-based answer, I'm not familiar with arrays but would love to learn

Thank you! This worked but we're only getting one contigs.fasta file when we run. Should we be getting a contigs.fasta per sequence pair or one for all the sequences we ran?

I am not familiar with spades but going by M__'s answer, you need something at the -o param as well. I'll edit my answer but you should test running 2 samples at the same time before running all of the scripts concurrently.

The job array is a formal approach to solving this via SLURM
What I am not clear about is the number of files SRR8599992_1.fastq.gz SRR8599992_2.fastq.gaz I originally thought the final number of files was around 100, I've guessed at 120
#SBATCH --job-name=pangenome_spade
#SBATCH --account=lilianasalvador
#SBATCH --ntasks=120
#SBATCH --cpus-per-task=16
#SBATCH --mem-per-cpu=5gb
#SBATCH --time=02:00:00
#SBATCH --array=599992-5600100

module load python/3.9
module load spades/3.15.0-hpc

spades.py -1 'SRR'${SLURM_ARRAY_TASK_ID}'_1.fastq.gz' -2 'SRR'${$SLURM_ARRAY_TASK_ID}'_2.fastq.gz' -o /groups/lilianasalvador/noah_100/'SPAdes_output_'${SLURM_ARRAY_TASK_ID}'.fastq' -t 16 --careful

I would test this via an echo test script writing to a specific file to make sure the spades.py command produces the correct command line for each iteration.
Things to note ...

Array jobs is specified by the $SLURM_ARRAY_TASK_ID variable and #SBATCH --array=599992-5600100
ntasks is the total number of spades assemblies you are performing here I guessed at 120
You are parallelising via -t 16 therefore this must be reflected in your submission job via --cpus-per-task
The ouput file MUST use $SLURM_ARRAY_TASK_ID otherwise you will continually override your out file.

This is not an optimal SLURM job because ...

-t 16 you need to check with your HPC engineer, this might take a long time to get through the queue and an alternative architecture to permit 16 CPUs can get onto the cluster. HTC users will jump ahead of you, so on a hypervisor setup you can specify a separate cluster for your work. This is not well explained but hope it makes sense.
Shifting to high throughput computing (HTC) could be considered here to get this through a lot quicker, e.g. this is 120 sequential assemblies, 10 sequential assemblies for 12 separate submissions would be a lot quicker (you just alter the "array job" inputs range.