text:
Can't load files into PopART
I am trying to create haplotype networks in PopART (Population genetics with Reticulated Trees*) which I have successfully done many times in the past.
Two weeks ago the program was working fine for me but recently I am unable to load files. What happens is that when I go to load a nexus file through the "nex" button, or through the "file/import/alignment" or "file/open" options nothing loads.
There is no error message, instead nothing happens. This occurs with both files I have recently created and with files that I have successfully loaded in the past. I've sent files to other people and everyone has been able to load the files using PopART on their computers. I am running macOS 10.15.6 and have tried turning my computer off and on as well as uninstalling and reinstalling PopART. I also thought that maybe Java wasn't up to date on my computer but I don't think PopArt runs with Java. Any help would be very appreciated!

PopArt builds reticulated HAplotype networks such shown below

Answer from @user438383, converted from comment:
Hi, this seems like quite a niche problem that you are facing, perhaps a bug in the GUI of some kind. Although it's on-topic, I think you're unlikely to find somewhere here who can help (unless the popart devs lurk here). I would suggest emailing the authors, but they point towards the github which mostly seems to be ignored.

Thanks, ill email the authors.

PopArt have a mailing list popartusers+subscribe@googlegroups.com
If you go here you can subscribe to it. This sounds like an issue with Mac compatibility. OS 10.15.6 isn't a version I understand and suspect there's a typo. Currently I believe it is OSX 12.1, thus OS10.12.1
PopArt was written for OSX10.6 and upwards. HOWEVER, there can be issues with major system upgrades so Monterey OSX12, I think was a major revision. PopArt is no longer supported so it would be sensitive to changes in Mac OS.
My suggestion is to try this on a Windows machine, if you get the same error the suspicion is there has been a small change to the input format.


text:
Converting gene symbols between species
I have a list of genes in a mouse model
How I can find the equivalent gene symbol/name for human from them?
This is my list
CD45
MHC II
CD11b
Ly6C
Ly6G
F4/80
CD11c
CD38
Arg1
SiglecF
CD206
CD62L
CD103
iNOS
PD-L1
TNFa
CD64
TCRgd
Foxp3
RORgt
CD8α
Tbet
CD25
IFN-γ
CD44
CD86
CD80
PD-1
B220
NK1.1
CD19
CD4
TCR β

From these I only find two homologous

Easiest way is BioMart. Help video to get you started here. Use the mouse genes dataset and filter by the list of gene names, get the mouse gene name and the human gene name (listed under homologues) as attributes.

Thank you so much I checked and likely from 33 genes, two genes have corresponding please look at the screen shot in my main post? How I find what is the conventional (standard names of these genes please?

There is no "conventional name", you're going to have to be more specific. The table has both EnsEMBL Gene Identifiers and HGNC Gene Symbols. That should suffice for any analysis.

That's just the first ten lines of your output, use the View All or downloads to get the full table, as was stated in the help video I sent you. You may also wish to untick some of the default options, as you're getting a line in your table for every transcript, which you don't need.

source: Stack Exchange

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

sort: votes


text:
how to get normcounts for singlecellexperiment object?
I need to perform differential expression analysis using the scDD package from R, but I am not able to since I miss the normcounts assay in my SCE object (of course in the example they show, the assay is already there).
As I found out in the preprocessing steps, one can use this function to generate the missing assay, except that it does not work:
# this is the SCE object
> cTSubSet
class: SingleCellExperiment
dim: 19767 149
metadata(2): Samples Samples
assays(2): counts logcounts
rownames(19767): Xkr4 Gm19938 ... CAAA01147332.1 AC149090.1
rowData names(4): ID Symbol Type chrLoc
colnames(149): AAAGGATGTTCACGAT-1 AACCATGTCTTTCGAT-1 ...
  TGTGGCGTCTCTTAAC-1 TTTGGTTAGACGCTCC-1
colData names(17): Sample Barcode ... slm_1 slm_1.5
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

# setting parameters as suggested in the vignette linked above
> paR <- list(alpha=0.01, mu0=0, s0=0.01, a0=0.01, b0=0.01)
# setting zero.thresh to zero since I already removed unwanted genes
> cTSubSet <- preprocess(cTSubSet, zero.thresh=0, scran_norm=TRUE)

# the not so useful error:
Performing scran Normalization
Error in t.default(Data) : argument is not a matrix

I was also looking for other possible ways of obtaining such normcounts assay and I found this description:

normcounts: Normalized values on the same scale as the original
counts. For example, counts divided by cell-specific size factors that
are centred at unity.

Still, I can't figure out what this may mean. Is there a way of obtaining such normcounts?

You might as well try monocle3::normalized_count.
It takes SingleCellExperiment object.

yep that is a possibility. I actually found another one using edgeR, i.e.
y <- DGEList(counts=(counts(cTSubSet)), group=cTSubSet$condition)
y <- calcNormFactors(y)
however, I am still working on this

I do not know this package but it seems to call scran so why bothering and not call scran directly, see for an extended discussion + code:
https://bioconductor.org/books/release/OSCA/normalization.html#normalization-by-deconvolution

yes, I am familiar with the OSCA pipeline. At this point I am just wondering if the values that are returned by computeSumFactors() can be considered as normcounts.

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


text:
Programmatically adding hydrogen and remove water to multiple PDB files
I'm trying to dock a ligand to several hundred PDB files (receptors).
I thus need to prepare the those PDB files like removing water and adding hydrogen.
I can do that manually using Autodock or Pymol. But it's too troublesome.
Is there a way to do it programmatically with Python or R or any command line procedure?

If you did it manually in PyMOL, it can also be scripted. But note that adding hydrogens is not that simple. There is more than one way of deciding which hydrogens are present (you may consider pH) and what are their coordinates. So I'd first think what method of adding hydrogens to use.

@marcin In pymol I use h_add command. That is not enough? How would you go about scripting it?

https://chemistry.stackexchange.com/questions/129157/adding-hydrogens-to-multiple-pdb-files-using-openbabel-gui-software here they use openbabel, not scripting but in batch

I can't tell if that's enough for what you're doing. h_add in Python script would become pymol.cmd.h_add('all')

Michael Word's reduce is an established method to do that

GROMACS can add hydrogens. And  water you can filter out with cat filename | grep -v ...

I believe PDB2PQR CLI will do the work wonderfully. Don't let the name trick you: PQR files are organized like PDB ones. Under the hood it runs propka, which is state-of-the-art for predicting a protein residues protonation state. The best part is that you can use PDB2PQR web server and they will give you the corresponding CLI arguments for each option you choose. It also automatically removes water (and other non-protein molecules).
You can also try looking into tutorials on how to setup proteins for dynamics on GROMACS, AMBER or OpenMM and stumble on their own tools for protonating stuff (respectively pdb2gmx (CLI), protonate (CLI), and modeller.addHydrogen (python API)), but as far as I know propka computes residues pKi in a knowledge-based manner and out-performs physics-based (if not just totally automatic, without guess) approaches like these.
For some context: acidic or basic amino-acids in proteins may have different pKi depending on their local environment (i.e. if they are in contact with solvent or not, and residues surrounding them), which is why "adding hydrogen to proteins crystallographic structures at a constant pH" is not as easy as it sounds.


text:
Illumina vs PacBio vs Nanopore price
Can someone help me with comparing the cost per run of Illumina, PacBio and Nanopore sequencing technologies?

Voted to keep open. I honestly think, although it's a wide question, it's essential for planning of any genomics project and it's not like these comparisons don't exist!

Could you please add more context to this, or narrow your question? Searching Google Scholar for illumina nanopore pacbio turns up thousands of comparative results papers, and StackExchange doesn't work well with thousands of answers.

You don't think that the kind of experiment you want to run is at all relevant?

voted to keep open, agree with Kamil S Jaron. This is information that is often intentionally obscured to the consumer by sequencing providers, but is essential for anyone trying to plan. compiling information is important for the community.


text:
Getting conventional gene symbol for Seurat
I have a Seurat object made by human single cells
I am mapping some genes on that but no sign of expression
When I GOOGLE for those genes I see the genes have different names
How I know Seurat uses which gene name please
This is my list
CD45
MHC II
CD11b
Ly6C
Ly6G
F4/80
CD11c
CD38
Arg1
SiglecF
CD206
CD62L
CD103
iNOS
PD-L1
TNFa
CD64
TCRgd
Foxp3
RORgt
CD8α
Tbet
CD25
IFN-γ
CD44
CD86
CD80
PD-1
B220
NK1.1
CD19
CD4
TCR β

From this list only 6 genes being mapped on the Seurat data

You have already asked this question and it has been answered - with the exact same response. Why are you opening multiple threads instead of exploring more on your own?

You might want to use Biomart to figure out the official names for those genes based on their aliases.

Can I trust NCBI? For instance this link which says TBX21 is an official gene name? https://www.ncbi.nlm.nih.gov/gene/30009

Yes, you can trust NCBI. Like I said in an earlier comment, HGNC provides the "official" gene symbols, and NCBI/EnsEMBL cross-referencing is reliable.


text:
querying gdc metadata issues with resulting fields
I have been trying to access metadata from the Genomic Data Commons; however, not all of the information I want is showing up in the resulting tsv file. How do I get sample type and analyte type to show up in the tsv file? I have tried altering what I pass as a field (towards the top of the code), or adding a different filter, but nothing seems to change what shows up in the tsv file. Any help is appreciated!
import requests
import json

cases_endpt = "https://api.gdc.cancer.gov/cases"

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = [
    "sample_ids",
    "submitter_sample_ids",
    "files.cases.sample_ids",
    "files.cases.submitter_sample_ids",
    "files.cases.samples.sample_id",
    "files.cases.samples.sample_type",
    "files.cases.samples.sample_type_id",
    "files.cases.samples.submitter_id",
    "samples.sample_id",
    "samples.sample_type",
    "samples.sample_type_id",
    "tissue_source_site.bcr_id",
    "tissue_source_site.code",
    "tissue_source_site.name",
    "tissue_source_site.project",
    "tissue_source_site.tissue_source_site_id",
    "diagnoses.classification_of_tumor",
    "diagnoses.tumor_stage",
    "diagnoses.tumor_grade",
    "files.cases.diagnoses.tumor_grade",
    "files.cases.diagnoses.tumor_stage",
    "files.cases.diagnoses.classification_of_tumor",
    "files.cases.samples.tumor_code",
    "files.cases.samples.tumor_code_id",
    "files.cases.samples.tumor_descriptor",
    "files.cases.samples.portions.slides.percent_tumor_cells",
    "files.cases.samples.portions.slides.percent_tumor_nuclei",
    "samples.tumor_code",
    "samples.tumor_code_id",
    "samples.tumor_descriptor",
    "samples.portions.slides.percent_tumor_cells",
    "samples.portions.slides.percent_tumor_nuclei",
    "files.cases.samples.oct_embedded",
    "samples.oct_embedded",
    "files.cases.samples.intermediate_dimension",
    "files.cases.samples.longest_dimension",
    "files.cases.samples.shortest_dimension",
    "samples.intermediate_dimension",
    "samples.longest_dimension",
    "samples.shortest_dimension",
    "files.cases.samples.portions.is_ffpe",
    "samples.is_ffpe",
    "samples.portions.is_ffpe",
    "files.cases.samples.pathology_report_uuid",
    "samples.pathology_report_uuid",
    "exposures.weight",
    "files.cases.exposures.weight",
    "files.cases.samples.current_weight",
    "files.cases.samples.initial_weight",
    "files.cases.samples.portions.weight",
    "samples.current_weight",
    "samples.initial_weight",
    "samples.portions.weight",
    "files.cases.samples.composition",
    "samples.composition",
    "files.cases.samples.freezing_method",
    "files.cases.samples.time_between_clamping_and_freezing",
    "files.cases.samples.time_between_excision_and_freezing",
    "samples.freezing_method",
    "samples.time_between_clamping_and_freezing",
    "samples.time_between_excision_and_freezing",
    "files.cases.samples.days_to_sample_procurement",
    "samples.days_to_sample_procurement",
    "files.cases.samples.preservation_method",
    "samples.preservation_method",
    "samples.days_to_collection",
    "files.cases.samples.days_to_collection",
    "files.cases.samples.portions.analytes.analyte_type",
    "files.cases.samples.portions.analytes.analyte_type_id",
    "files.cases.samples.portions.analytes.aliquots.analyte_type",
    "files.cases.samples.portions.analytes.aliquots.analyte_type_id",
    "samples.portions.analytes.analyte_type",
    "samples.portions.analytes.analyte_type_id",
    "samples.portions.analytes.aliquots.analyte_type",
    "samples.portions.analytes.aliquots.analyte_type_id",
    ]

fields = ','.join(fields)

filters = {
    "op": "in",
    "content": {"field": "cases.project.project_id",
                "value": ["TCGA-BRCA"]},
}

params = {"filters": json.dumps(filters),
          "format": "TSV",
          "size": "12000"}

response = requests.get(cases_endpt, params=params)
content = str(response.content, "utf-8")

with open("addingsampletype5.tsv", "w") as file:
    file.write(content)

When I visit the cases_endpt page you linked directly on Firefox v88, I get a filterable JSON file. Looking for "type" only comes up with "disease_type", not sample or analyte:
data
 hits
  0
   disease_type "Adenomas and Adenocarcinomas"
  1
   disease_type "Ductal and Lobular Neoplasms"
  2
   disease_type "Ductal and Lobular Neoplasms"
  3
   disease_type "Adenomas and Adenocarcinomas"
  4
   disease_type "Squamous Cell Neoplasms"
  5
   disease_type "Adenomas and Adenocarcinomas"
  6
   disease_type "Squamous Cell Neoplasms"
  7
   disease_type "Adenomas and Adenocarcinomas"
  8
   disease_type "Squamous Cell Neoplasms"
  9
   disease_type "Ductal and Lobular Neoplasms"


text:
alternative to HOMER's mergePeaks
Recently, our lab ran a 192 sample experiment through our ATACseq pipeline. In doing so, HOMER's mergePeaks told us that our 512 GB RAM server had too little memory to process all of the samples at once. So here we are trying to replace it with an alternative.
It seems that bedTools merge would be the way to go to replace it. However, going into this step, each sample has its own bed file of peaks called by MACS2. The mergeBed only takes in 1 input file that must be sorted by chromosome and peak start position before we run it. I can run a cat call on all of the samples followed by a sort, but those individual sample bed files are already sorted. Is there a way to run mergeBed with multiple input files at once? I tried mergeBed -i 1.bed -i 2.bed > out.bed, but that is equivalent to mergeBed -i 1.bed > out.bed.
In image format, I want the leftmost example, where the three peaks come from three separate files instead of all being described in the same file:
Additionally, bedtools does have a poorly documented function called multiIntersectBedthat appears to run the cat and sort together, but is both slower and doesn't merge bookended sections--eg files look like this:
chr1 4658186 4658322
chr1 4658322 4658775
chr1 4658775 4658777
chr1 4671049 4671351

When they should look like this:
chr1 4658186 4658777 # This was 3 lines before
chr1 4671049 4671351

In summary, I have sorted bed files 1.bed and 2.bed from the same experiment, and I want to merge all bookended and overlapping peak sections between them. How do I do so in 1 step that doesn't involve HOMER's mergePeaks? Thank you!
EDIT: Here's the speed of my current solutions at n=192:
HOMER: Unstable

`multiIntersectBed | mergeBed`

real    3m53.230s
user    6m28.166s
sys     0m14.912s

`cat | sort | mergeBed`

real    0m24.860s
user    0m30.939s
sys     0m2.946s

`cat > tmp1; sort > tmp2; mergeBed > out; rm tmp?`

real    0m11.824s
user    0m20.710s
sys     0m2.021s

`sort -m | mergeBed`

real    0m17.905s
user    0m25.112s
sys     0m1.816s

`sort -m --batch-size=192 | mergeBed`

real    0m14.687s
user    0m25.589s
sys     0m0.747s

`sort -m > tmp ; mergeBed -i tmp ; rm tmp`

real    0m17.472s
user    0m15.251s
sys     0m2.212s

So... Followup questions, why does sort -m take longer than sort | cat especially when the sys time is indeed smaller? And why does piping take more time than creating temporary files?

Via BEDOPS, you could combine bedops --merge with process substitutions sorting N files:
$ bedops --merge <(sort-bed 1.bed) <(sort-bed 2.bed) ... <(sort-bed N.bed) > answer.bed

The peak memory usage of a typical bedops --merge run should be around 10MB, if I remember correctly, as a result of using sorted inputs.
The memory usage of each sort-bed process can take up the size of each file. You can specify --max-mem and --tmpdir options to instead use disk space to sort files. You may want to do this if the sum of file sizes is equal to or greater than ~512GB, or however much system memory is available (typically less).
If your inputs are already sorted, or if your ATACseq pipeline is set up to create sorted files, then sort-bed does not need to be run again, and you can just use the merge operation on its own:
$ bedops --merge 1.bed 2.bed ... N.bed > answer.bed

The memory usage is then only what is needed for merging, as discussed above.
Using bedops --merge will merge adjacent ("bookended") intervals, so long as they are in sort-order, e.g.:
$ echo -e "chr1\t4658186\t4658322\nchr1\t4658322\t4658775" | bedops --merge -
chr1    4658186 4658775

Notes:

The file answer.bed will be sorted per sort-bed for downstream operations, such as bedmap association of merged regions to peaks. Additional sorting is not needed.
Using sort-bed will often be faster than Unix sort on BED files.
The limit of N on most platforms will be 1021 files. Most OSes have a limit of 1024 file descriptors, and three are reserved for standard IO. This soft limit can be dialed upwards on some hosts where the user has admin access.

References:

https://bedops.readthedocs.io/en/latest/content/summary.html#bedops
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#merge-m-merge
https://bedops.readthedocs.io/en/latest/content/reference/file-management/sorting/sort-bed.html


text:
How to make 3d model of a protein that not exist in PDB?
I'm working on nsSNPs(Non-synonymous single nucleotide polymorphisms) of protein AKAP13 whose PDB not available.
I did one mutation in this protein and I want to show the effect by 3d modeling without using a PDB but don't know how to do.
Can anyone help Please!!

Uh-oh. You seem to have been downvoted twice because your question was a how do I go about do X?, yet your reply asking for more info showed you already had tried things: the correct question should have been a I tried A & B, but cannot do X in the first instance.

Note to passersby: do not downvote as the user is a newcomer.

No PDB
Getting a crystal structure is hard work for crystallographers, even with a high throughput systems, so they make constructs with only known domains or regions of particular interest. What happened here is that only two parts have been cloned, expressed and crystallised. The rest may be crystallisable or solvable by cryoEM, but may be too disordered.
Homology modelling
If the protein has a close homologue it can be threaded, that is it's sequences is "threaded" over the structure of the homologue. Two site provide these ready-made.

SwissModel
ModBase

The former is very clear, has an NGL viewer, a handy circular map and the assessment is clear. For your protein it is showing four existing crystal structures, but these cover only a small part of it. The homology modelled section is small and not good here, but handy as the model has two numbers Qmean of –4 and identity of 23%, two red flags. So the next place to check are the following servers. They thread against different structures and hybridised. Missing parts are predicted ab initio (from the scratch). The latter parts look like spaghetti and cannot be trusted.

Phyre2. <1 day turnover. Do pay attention to the colour rails. If your region of interest in low confidence and is based on nothing, then do not trust it at all. Also there it is good to double check keep an eye on the sequence alignments of your position as misalignments can happen.
I-Tasser. 3-5 day turnover. The CAPS server winner year after year. You 5 models representing each cluster. If they look very different and the metric is bad for the top, then you have very low confidence and must look at the all if even.
EVFold. The amino acid co-ocurrence is rather interesting
AlphaFold... Google keeps it secret. Nobody except them can use it.

Here is an example of a spaghetti model (not your protein):

The except for the red bar, which is modelled against homologues and forms a nice domain within the spaghetti ball at the bottom, all of the rest is ab initio and should not be trusted as the warning in the "Multi-template and ab initio information" section says.
A nice addendum to a model is Consurf, which maps conservation to the models as B-factors. In a crystal structure, the core is conserved and the surface is not. Metrics are nice, but human insight is better.
A big problem are protein that are part of a complex, because alone they are not structured and their structure is dependant on a neighbour. Models are great if they are threaded against a model that forms the same complex, eg. mouse to human, but not if they are hybridised.
In the case of multidomain protein that are large. It is best to split the protein up. Furtehrmore, I-Tasser and Phyre have length limits because the longer the protein to predict is, the higher the chance of error.
One can determine where the domain boundaries may be —a Blast or a few iterations of Psi-blast will show conservation boundaries nicely. PFam has a lot to say about the protein architecture: enter link description here
Regions of high disorder cannot be and should not be modelled as they do not form a structure free in solution —unless wrapped around another protein as mentioned.
Wisdom to step away
However, the real question is not how to but should one. If one only get ab initio results, do not trust them at all —count your losses as opposed to embarking on a castle in the sky. For a SNV, check if it is close to a post-translational modification in PhosphoSitePlus or if there is a linear motif in ELM.

Answer from @o-laprevote, converted from comment:
Your protein includes intrinsically disordered regions (https://www.uniprot.org/uniprot/Q12802#family_and_domains). Showing a predicted structure may give a wrong idea of what it looks like: I'd consider simply showing the ordered part (of which a structure already exists) plus a "spaghetti" starting from its N-terminal residue.
In this case the ordered regions would be in the "interaction with ESR1" domain (1968-2338), which as you can see has been solved experimentally (PDB 4D0N), and the sequence spanning from 965-1431 (part of it has been solved, PDB 2DRN). What I mean by "spaghetti" is basically: render these two regions then draw curvy arabesques/lines between, before and after them to show it is disordered (starting from their N-terminal or C-terminal).
Such a cartoon is best done in a program that is shows it utterly differently —to not confuse the viewer as this is for display purposes only. Namely, render the parts in PyMOL, but then add the connecting lines in Photoshop or Powerpoint or MS-Paint/Preview.

thank you for your detailed explanation, but the region that I have to visualize is 2665 – 2813 and it is disordered, and can't view it in Pymol because I don't have PDB for this region


text:
Downloading PSSM file from NCBI PSI-Blast results after the UI update
It used to be easy to download off a results page of a PSI-Blast from NCBI the PSSM (scorematrix). Now I cannot seem to find an asn file that is a scorematrix —Download>ASN.1 is the output of the search. Whereas before it was in "format results" IIRC.
I was after this to check something in the PSSM viewer or do further cycles of PSI-Blast with different settings —the latter especially handy given that at present the exclude taxa checkbox is troublesome. Where is the option gone?

Addendum: a blast search with only excluded taxa never used to glitched but now does, but adding biota taxa prevents it from glitching.

Have you tried clicking on the "Back to Traditional Results Page" link in the top right corner? If you're expecting a different interface, that might help work out where to find the right information.


text:
Best way to dock subset of protein PDB and definition of a ligand
I have extracted residue 45 to 88 from a PDB file (2YRQ).
With that I have a new smaller PDB file. Let's call it 45_88.pdb
And I want to dock it to some other full protein.
My question is what's the best way to dock it. I'm considering Autodock Vina.
But I'm not sure if it's appropriate since it is for ligand-protein docking.
What's the definition of a ligand? Can a protein with only 42 residues be considered a ligand?

I would advise using the term peptide more so than ligand, which also includes, and generally is used to refer to, small molecules.

Answer from @pippo1980, converted from comment:
https://pubs.acs.org/doi/10.1021/acs.jctc.9b01208 Comprehensive Evaluation of Fourteen Docking Programs on Protein–Peptide Complexes


text:
3d model of a protein via nucleotide sequence
I really appreciate @Matteo Ferla's detailed answer to my previous question.
I've already tested several tools to know if a mutation can cause a modification in the function of protein or not. There are some tools such as SIFT, PolyPhen 2, PROVEAN, SNPs&GO, I-Mutant, and so on. These tools give us the possibility to check if this mutation is deleterious or tolerant, So after using my own data on these tools I found that most of the tools give me ‘deleterious ‘ meaning that the structure of my mutated protein will be changed. So, my problem is about the mutated protein that does not have a PDB Code which I can use in Pymol, mcSm, SDM, DUET…. to visualize the 3d structure I tried to use BioPython but it requires a PDB code.
You talked about the swissmodel, I had tried to use this server but it gives me the same PDB code for the original protein and the mutated protein with a similarity of 100% and QMEAN approximately the same, as shown in the following figures:

After that, I had used this PDB code in PyMol but I didn't find the position 2700 in which the mutation occurs.
The ModeBase displays the same PDB code as Swiss-Model in which there is no 2700 position.

The I-TASSER and EVFold require a small amount of data, since I am working on position 2700 these tools are not useful for me.

There is no method to visualize a 3d model with just a FASTA file or how to do my own PDB file based on a FASTA file using Python or R or any other language?
I have Fasta format data (not PDB code) which I would like to visualize in 3d structure using R or Python or anything that can help me because I can't use tools of 3d structure(i have not PDB code).
I've already tested several tools to know if a mutation can cause a modification in the function of protein or not. There are some tools such as SIFT, PolyPhen 2, PROVEAN, SNPs&GO, I-Mutant, and so on. These tools give us the possibility to check if this mutation is deleterious or tolerant, So after using my own data on these tools I found that most of the tools give me ‘deleterious ‘ meaning that the structure of my mutated protein will be changed. So, my problem is about the mutated protein that does not have a PDB Code which I can use in Pymol, mcSm, SDM, DUET…. to visualize the 3d structure I tried to use BioPython but it requires a PDB code.
As long as the mutation is not known yet, I couldn't use these tools because they need the PDB code that I haven't.
I would be grateful to those who will help me!!


text:
Connection between Detected Genes and The Read Counts
I have been trying to understand the Seurat for analysing scRNA-seq data. It comes to my mind that the main data is organised in the Seurat object with rows as genes and columns as the cells, and the entries are the gene read count.
I am confused with the statement, "the connection between the detected genes in the cells and the read counts".
What is exactly the meaning of the detected genes in the cells biologically and its connection to the read counts? Can we say if the read count of gene A is 10 it means this gene is detected 10 times in a cell? How to interpret the detected genes biologically?

This needs more focus and a reference / context of that quote.

Answer from @ATpoint, converted from comment:
A detected gene is one that has a count > 0, and generally the deeper you sequence the more genes you detect until you reach a saturation which is being influenced by the number of biologically-expressed genes and the single-cell technology you're using, e.g. droplet vs plate-based and the library setup, e.g. end-tagged vs full length.
Counts cannot influence the nature of gene translation. A count is an observation; a translated gene sequence is a physical thing. Whether or not a gene can be translated into protein has no relationship to how often it is observed in existence.
Even with count zero (i.e. no detection in scRNA-seq), it does not mean that the gene is not present in the cell. NGS technology is a relative quantification, I doubt it is meaningful to try to directly make statements about the connection between read counts and biological meaning in the cell.


text:
Count of each sequence length from a fasta file with header using len function of python
I have a fasta file with 18 sequences inside. I want to extract the length of each sequence with header using the len function. My file looks like:
>DR-1
ATGTATGAGGGCAAGAAGACGAAGAACATGTTCCTGACCCGGGCCCTGGAGAAGATTCTG
GCCGACAAGGAAGTGAAGAAGGCGCATCACTCCCAGCTGCGCAAAGCGTGCGAGGTGGCC
CTAGAGGAAATAAAAGTTGAAACTGAGAAACAAAGCCCTCCTCATGGAGAAGCAAAAGCT
GGATCAGGCACCCTTCCACCAGTGAAATCAAAGACAAATTTTATTGAAGCAGACAAGTAC
TTTTTACCCTTTGAGTTGGCATGCCAGTCCAAATGTCCCCGCATAGTTAGTACATCTCTA
GATTGTTTGCAGAAACTTATTGCTTATGGGCACTTGACTGGCAGGGCTCCAGATAGTACG
ACACCAGGCAAAAAATTAATTGATAGAATTATTGAAACAATATGTGGCTGCTTCCAAGGT
CCTCAAACAGATGAAGGAGTTCAGTTGCAGATAATAAAGGCTTTACTTACCGCAGTGACA
TCACAACACATAGAAATCCATGAAGGAACTGTGCTGCAAGCTGTGAGGACATGCTACAAT
ATCTATTTAGCCAGCAAAAATCTCATCAATCAGACAACAGCCAAAGCTACTCTTACTCAA
ATGCTGAACGTTATCTTTGCACGCATGGAAAACCAAGCATTACAGGAGGCCAAACAAATG
GAAAGAGAAAGACATCGGCAGCAGCAGCATTTGTTACAGTCTCCAGTAAGCCATCATGAG
CCCGAGTCACCTCATCTTAGATATTTGCCACCTCAGACTGTTGATCATATAAACCAAGAG
CATGAAGGGGATCTTGAGCCCCAAACACATGATGTGGATAAAAGTCTCCAAGATGATACA
GAGCCTGAAAACGGATCTGATATTTCCAGTGCAGAAAATGAACAGACAGAAGCTGATCAG
GCAACTGCAGCTGAAACATTATCTAAAAATGATATCTTGTATGATGGAGACTATGAAGAA
AAGCCACTTGATATTGTCCAGAGCATTGTAGAAGAAATGGTGAACATCATTGTTGGAGAC
ATGGGAGAAGGGATGGCTATCAGTGCAAGTACAGAAGGCAACACTGGAACTGTAGAAGAT
GGTAGTGACAGTGAAAACATCCAAGCAAATGGAATCCCAGGAACACCCATTTCTGTTGCA
TATACACCATCGTTACCTGACGACAGGTTATCTGTCTCCTCCAATGATACTCAGGAATCT
GGAAATTCTTCAGGACCTTCACCTGGTGCTAAGTTTTCCCACATTTTACAAAAGGATGCC
>SK-1
TTTCTAGTGTTCAGGTCATTGTGTAAGCTGTCAATGAAACCACTGTCAGATGGACCACCA
GATCCAAAGTCTCATGAGTTACGGTCCAAGATTCTTTCATTGCAGTTGCTTCTGTCCATT
CTGCAGAATGCAGGACCTGTTTTCAGGACAAATGAGATGTTTATTAATGCTATCAAGCAG
TATCTATGTGTTGCACTCTCAAAAAACGGAGTCTCATCTGTTCCAGAGGTTTTTGAGCTT
TCTCTTTCCATATTTCTTACTTTGTTGTCAAACTTCAAGACACATCTGAAGATGCAAATT
GAGGTATTCTTTAAAGAAATCTTCCTGTACATTTTGGAAACGTCTACAAGCTCATTTGAC
CATAAGTGGATGGTTATTCAGACACTGACAAGGATCTGTGCAGATGCTCAGAGTGTGGTG
GACATTTATGTAAATTATGACTGTGACCTAAATGCAGCAAATATATTTGAAAGACTAGTA
AATGACCTGTCGAAAATTGCCCAAGGAAGGGGCAGCCAAGAACTTGGTATGAGTAATGTT
CAGGAATTGAGCCTGCGAAAAAAGGGGTTAGAATGTTTAGTTTCAATTTTGAAATGTATG
GTTGAATGGAGCAAGGATCAGTATGTGAACCCCAACTCTCAAACTACCCTTGGTCAGGAA
AAACCCTCAGAGCAAGAGATAAGTGAAGTTAAACACCCAGAGACTATAAACAGATATGGA
AGTTTAAACTCCTTAGAATCCACATCATCATCAGGAATTGGCAGCTACAGTACACAGATG
TCTGGCACTGATAATCCAGAACAATTTGAGGTCCTAAAGCAGCAAAAAGAAATAATAGAA
CAAGGGATTGATTTATTTAATAAGAAACCAAAGAGAGGAATACAATACCTCCAAGAACAA
GGGATGCTTGGCACCACACCCGAAGATATAGCTCAGTTCTTACATCAAGAAGAAAGATTA
GACTCTACTCAGGTGGGTGAGTTCCTGGGAGATAATGATAAATTCAATAAAGAAGTCATG
>SC-1
TATGCATATGTAGACCAACATGACTTTTCAGGAAAAGACTTTGTTTCAGCTCTTCGTCTA
TTTCTGGAAGGATTCCGGCTTCCAGGAGAAGCTCAGAAAATTGATCGATTAATGGAAAAA
TTTGCTGCAAGATACCTAGAGTGCAACCAAGGACAAACCCTTTTTGCTAGTGCAGATACT
GCTTATGTTTTGGCCTACTCAATTATCATGTTGACTACAGACCTTCACAGTCCACAGGTT
AAAAACAAAATGACAAAAGAACAGTACATTAAAATGAATCGTGGTATCAATGACAGTAAA
GATCTCCCTGAAGAGTATCTATCAGCCATTTATAACGAAATAGCTGGGAAAAAAATATCA
ATGAAAGAAACAAAAGAACTAACAATCCCTACAAAGTCAACTAAACAGAATGTAGCCAGT
GAAAAGCAAAGACGACTTCTGTATAACTTAGAAATGGAACAGATGGCCAAGACGGCTAAA
GCACTCATGGAAGCTGTGAGCCATGTTCAGGCACCTTTTACAAGTGCAACACACTTGGAG
CACGTGAGGCCGATGTTTAAGTTGGCTTGGACACCTTTTCTGGCTGCATTCAGTGTGGGT
CTCCAAGACTGTGATGATACTGAAGTAGCCTCTCTTTGCCTAGAAGGGATAAGATGTGCT
ATCAGGATTGCATGCATTTTCAGCATTCAGCTGGAAAGAGATGCATATGTCCAAGCACTA
GCCAGATTTACCTTACTTACAGTGAGTTCTGGTATTACTGAAATGAAACAGAAGAATATC
GACACAATAAAAACACTCATCACGGTGGCTCACACTGATGGAAATTATTTAGGAAATTCG
TGGCATGAGATCCTGAAGTGCATCAGTCAGTTAGAGCTGGCACAACTTATAGGAACTGGA
GTGAAGCCTCGGTACATTTCTGGGACTGTGAGAGGCAGAGAAGGGTCTCTCACTGGGACA
AAAGATCAGGCTCCAGATGAATTTGTTGGTCTAGGGCTAGTTGGAGGAAATGTGGACTGG
AAGCAGATAGCCAGTATTCAGGAGTCTATTGGAGAAACCAGCTCTCAAAGTGTGGTAGTT
GCTGTGGATAGAATATTCACCGGATCCACAAGGCTAGATGGAAATGCTATTGTTGATTTT
GTCCGTTGGCTATGTGCTGTGTCTATGGATGAATTGCTTTCTACTACACACCCAAGAATG
TTTAGTCTACAAAAAATAGTAGAAATATCATACTATAACATGGGAAGAATAAGACTGCAG
TGGTCTCGGATTTGGGAAGTTATTGGAGATCATTTTAATAAGGTTGGTTGCAATCCTAAT
GAAGATGTAGCCATATTTGCAGTAGACTCTCTGAGGCAGCTATCAATGAAGTTTTTAGAG
>SC-2
AAAGGAGAGCTTGCAAACTTCAGATTCCAAAAGGATTTCTTACGACCTTTTGAACATATA
ATGAAACGAAACAGGTCTCCAACAATTCGAGATATGGTTGTACGGTGTATAGCACAAATG
GTTAATTCTCAGGCTGCAAATATTCGTTCAGGATGGAAGAACATTTTCTCAGTATTCCAT
CTAGCTGCATCAGACCAAGATGAAAGCATAGTAGAACTTGCATTTCAGACAACAGGGCAC
ATTGTCACTCTTGTATTTGAAAAACACTTTCCGGCCACCATTGATTCTTTCCAGGATGCA
GTGAAGTGTTTGTCTGAATTCGCATGTAATGCAGCTTTCCCAGACACAAGTATGGAAGCA
ATCCGACTTATTCGCCATTGTGCAAAATATGTGTCTGATAGACCGCAGGCTTTCAAGGAG
TACACAAGTGATGACATGAGCGTGGCACCTGAAGACAGGGTGTGGGTGAGAGGGTGGTTC
CCGATCCTCTTTGAGTTGTCTTGTATCATCAATAGATGCAAATTAGATGTAAGAACCAGG
GGCTTAACAGTAATGTTTGAAATAATGAAAACATATGGCCACACTTATGAAAAACATTGG
TGGCAGGACTTATTCAGAATTGTTTTCAGAATCTTTGACAACATGAAATTGCCAGAACAG
CAGACAGAGAAAGCTGAATGGATGACAACAACTTGTAACCATGCCCTTTATGCAATATGT
GATGTATTCACCCAGTATTTAGAAGTTCTCAGCGATGTACTTTTGGATGATATTTTTGCC
CAGCTTTACTGGTGTGTACAACAAGACAATGAGCAGTTAGCACGATCTGGTACAAACTGC
TTAGAAAACGTGGTTATTCTGAATGGTGAAAAGTTTACTCTTGAAATCTGGGATAAAACG
TGTAACTGCACCCTGGATATCTTCAAAACCACAATTCCTCACGCGCTTCTGACTTGGCGT
>TA-1
CCCACCTCTGGCGAGGCCGAACCCCCATCTCCATCAGCTGTGAGTGAAAAGCCGTTGGAT
GCAATATCACAAAAATCTGTAGACATACATGATTCTATTCAACCAAGATCTTCAGATAAC
AGACAGCAAGCACCACTGGTTTCTGTTTCTACTGTGAGTGAAGAAGTTAGCAAAGTGAAA
TCCACAGCAAAATTCCCAGAACAAAAATTGTTTGCTGCCCTGTTGATTAAATGTGTTGTG
CAGCTGGAACTCATCCAGACCATCGACAACATTGTCTTCTTCCCAGCCACAAGCAAGAAG
GAGGATGCAGAGAACTTAGCTGCAGCCCAGAGAGATGCAGTAGACTTCGATGTTCGAGTT
GATACTCAAGACCAAGGAATGTACCGCTTTTTAACATCACAACAACTGTTCAAGCTACTG
GACTGCTTATTAGAGTCGCATAGATTTGCAAAAGCATTTAATTCCAACAATGAACAGAGG
ACTGCTCTGTGGAAAGCAGGCTTCAAAGGCAAGTCCAAGCCCAATCTTCTAAAACAGGAG
ACAAGCAGCCTGGCGTGTGGGCTGCGCATTCTCTTCCGGATGTACATGGACGAGAGCAGA
GTCAGTGCCTGGGAAGAAGTCCAGCAGAGGCTTTTAAACGTCTGCAGGGAAGCGCTGAGC
TACTTCCTCACTCTAACATCCGAAAGTCACCGGGAAGCCTGGACTAACTTACTGTTACTG
TTCCTCACTAAAGTTCTGAAGATCAGTGACAGCAGGTTTAAAGCTCACGCATCATTCTAC
TACCCTCTCTTATGTGAAATTATGCAATTTGACTTGATTCCGGAACTGCGTGCTGTTCTT
AGAAGATTTTTCCTGCGAATTGGAATAGTATTTCAAATATCACAACCACCTGAACAGGAA
CTTGGGATAAACAGGCAATGA
>DR-2
ATGCAGCCACCACCTCAGGCAGTCCCGTCTGGGGTGGCTGGACCACCCCCAGCTGGGAAT
CCTCGGAGCATGTTCTGGGCCAACAGCCCTTACAGGAAGCCAGCCAATAATGCACCAGTG
GCTCCCATAACACGCCCATTGCAGCCAGTAACGGATCCGTTTGCTTTTAATAGACAGACT
CTTCAAAACACACCAGTAGGAAGTTCATCCAAAAGCAGCCTGCCCAATTTGCCTGGCCCA
GCCCTCTCTGTTTTTTCTCAGTGGCCTGGTTTGCCTGTGACTCCCACAAATGCTGGGGAT
AGCTCCACAGGATTGCATGAGCCTCTGTCAGGAACTCTGTCACAGCCTAGAGCAGATGCC
AGTCTATTTCCTCCTGCATCTACCCCTTCTTCACTCCCTGGGCTGGAGGTGAGCAGGAAT
GCTGAGGCTGATCCCAGCTCAGGACATGAAGTCCAGATGCTGCCACATTCAGCTCACTAC
ATTCCAGGAGTGGGTCCTGAGCAGCCTCTCGGGGGCCAAATGAATGACAGTGGGTCCGGG
CCTGACCAACCCATGAATAGGCATGCCCCCCATGATGGTGCCGTGACCCACGCAGCATCT
CCTTTCCTCCCTCAGCCACAGATGCCAGGTCAGTGGGGGCCTGCACAAGGAGGCCCTCAG
CCCTCCTATCAGCATCACTCACCCTACCTGGAGGGACCTGTTCAGAACATGGGGCTCCAG
GCTGCCAGCTTGCCCCATTTTCCCCCTCCATCTAGTCTACACCAGGGACCTGGCCATGAG
TCTCATGCCCCACAAACATTTACTCCAGCATCCTTGGCCAGTGGTGAAGGGAATGAGATA
GTTCACCAGCAAAGTAAAAATCACCCATTGAGTAGCTTCCCCCCTAAACACACTTTTGAG
CAGAATTCCAGAATTGGGAATATGTGGGCAAGCCCAGAGTTGAAGCAGAATCCGGGAGTG
AATAAAGAGCATCTGCTAGACCCTGCCCATGTCAATCCCTTCACTCAGGGAAATAGTCCT
GAAAACCAGGCGCATCATCCCCCAGTGGCAGCAACTAACCATGCCCTGCAAGAAGCAGCC
TCTGGAGCTCTCTCAATGTTTTTTCAAGGGGAAGAGACAGAGAATGAGGAAAATCTCTCA
TCTGAAAAAGCAGGCCTCGATAAGAGATTGAACTTGGATAGCTTTTCTTCTACCTCCAGA
CTGGGCCATCCACCACCACCTGGAGCTAGTGGTGTCTACCAGGCCTTTCCTAGAGGCCCT
AGCAGTGAGGCTGCGCAGGAAGGAGATGCGCAACCTTATTTTTCGCAGTCTGTGGGCGTC
CGTCTTGACAAACAGAGCACTGTGCCCCCTGCTAATGATGCATGGGGTGATGTTCCGGGG
ACCGGGACCCGCTGTGCCAGTGGGCCACAGTGTGAGAATGTTGAGAACTTAGAATTTGTT
CAGAATCAAGAAGTTTTGCCCCGTGAGACCCTGAGTGTGGATCCCTTCCCCCTGAGTGAT
CAGATCCGATATGGGCCCCTTCCTGGGCCAGCTGCCTCCAGGCCTGCTACTGTGGGCCTC
ACCAGAGGTGGGGGGCTAAATCTTGAGGCCCCAGACACACCACTGCACCCTACACGACCT
GACAGTGTGTCATCCAGCTACAGCAGCCATAGCCACAGAAGTCCTCCTGGGTCAGCCAGA
CCCCAAGAATTGGTGGGCACATTCATTCAGCAAGAAGTTGGGAAACTTGAAGATGATACT
TCAGGGAGTTTTTTTAAGCAAATTGATTCTTCTCCTGTTGGAGGTGAAACAGATGAGGTG
ACTGGGAGCCAAAATTGCTGCAGCAGTCTGTCCCAGCCTTCAACCCCAAGCCCTCCAAAA
CCTACAGGAGTATTCCAAACAAGTGCAAACAGTTCTTTTGAGCCAGTAAAATCCCACTTG
GTTGGAGTGAAACCAGTTGAAGCAGATCGTGCCAACATGGTGGTTGAAGTGAGAGGGACC
CAATACTGTCCGAAGAAGCGCAGGGCAGCTGTTGCACCACCTGATGCTACCTCAGGAAAC
CTGGAACAGCCACCAGACAACATGGAGACCCCGTGTGCACCTCAGGCCTGCCCTTTGCCT
CTTAGCACCACTGGAGAGGCTGGGCAACTGGTTTCAAACACGGCAGGGACACCCCTGGAT
ACTGTGCGCCCAGTGCCTGATAAGAGGCCCTCAGCCAGGGCCCAGGGACCTGTGAAGTGT
GAGAGCCCAGCAACAACTCTGTGGGCACAAAATGAGCTGCCAGATTTTGGTGGCAATGTC
CTTCTAGCCCCTGCTGCCCCTGCACTTTATGTTCCTGTAAAACCAAAACCTTCAGAGGTT
GTTCATCATCCAGAAAAGGGAATGTCTGGACAGAAGGCCTGGAAACAAGGGTCTGTGCCA
>TA-2
CCTTTGCAGAACCAAGACCCTCCTGGTGCTTCTGAGAACCTTGAGAACCCTCCCAAGGTA
GGAGAAGAGGAGGCCCTTCCAGTGCAGGCAAGTTCTGGTTATGCGAGTTTGCTTTCTTCG
CCACCCACTGAATCTTTGCACAATCAACCTGTCCTGATTGCTCAGCCTGATCAAAGCTAT
AATTTGGCTCAGCCTATTAATTTTTCTGTGTCCTTACTGAATCCTAATGAGAAGAATCAG
TCCTGGGGAGATGCTGTGGTGGGGGAGAGATCCATAGTGAGCAACAACTGGGCTCTCGGT
GGTGATCCTGAAGAACGTGCTGCTCTGTCTGGAGTTCCAGCCAGTGCTGTCACTGGCGCG
TCTTTGCCTAGCAGTATCCCCCAGAATTGTGCCCCTCAAGGTTCGGGGTCTTCTGAAATG
ATTGCTAGCCAGTCTGCCAGTTGGTTGGTTCAACAGCTATCCCCTCAAACCCCTCAAAGT
CCGCATCCAAATGCAGAAAAGGGCCCTTCCGAGTTTGTTAGTAGCCCTGCTGGCAATACA
AGTGTGATGCTAGTTCCACCTGCAAGCAGTACTTTGGTTCCTAACAGTAATAAGGCAAAG
CACTCTAGTAACCAGGAAGAAGCTGTTGGCGCCCTAGACTTCACATTAAATAGAACTTTG
GAAAATCCTGTGAGAATGTACAGCCCATCTCCTTCAGATGGCCCAGCTTCTCAGCAGCCT
CTTCCCAATCATCCCAGACAATCTGGGCCTGGGCTGCATAACCAGGATCATTTCTACCAG
CAGGTAACAAAAGATGCTCAGGACCAGCATAGACTAGAAAGAGCCCAGCCAGAGCTGGTG
CCTCCTCGGCCACAGAATTCTCCCCAAGTGCCCCAAGCATCGTGTCCAGAGCCTTCAAAT
CCTGAAAGCCCACCCACACAAGGACAGTCTGAAAGCTTGGCCCAGCCACCAGCAAGTCCA
GCTTCAGTCAACACAGGCCAGCTGCTGCCTCAGCCACCACAGGCTTCCAGTGCATCCGTT
ACATCTACCAACTCCAGCCAGGCAGCTGTGCGGTCAGAGCAGCTGTGGCTACACCCACCA
CCTCCTAATACTTTTGGCCCAGCACCTCAAGACTTGGCATCCTACTACTATTACAGACCC
CTGTATGATGCCTACCAGTCTCAGTACCCTTCACCGTACCCGTCAGATCCTGGGACAGCC
TCGCTGTACTACCAGGATATGTATGGCCTGTATGAGCCTCGTTACAGACCCTATGACAGC
TCAGCATCTGCTTATGCCGAGAACCACCGCTATTCTGAGCCTGAGCGGCCCAGCTCCCGT
GCAAGTCACTATTCCGACCAGCTTGCTCCCAGGCAAGGATATCCTGAAGGATACTATAAT
TCCAAGAGTGGATGGAGCAGTCACAGTGATTACTATGCAAATTATTACTCTGGCCAGTAT
GATTATGGAGACCCCAGCCGCTGGGATCGTTACTATGGGTCTCGCCTTAGGGATCCTCGC
ACCTGGGACCGGAGGTACTGGTATGACTCTGAACATGACCCATACAGGAAGGACCACTAT
GCTTACAGTGACAGGCCTGAGAAATGTGATGATCACTGGAGGTATGACCCTCGCTTCACT
GGGAGCTTCGACGATGACGCTGAGATCCACAGGGACCCCTATGGAGAAGAAGCAGACAGA
CGCAGCATCCACAGTGAGCACTCGGCACGGAGCCTGCGCAGCACTCACAGTCTGCCCAGC
CGCCGCAGCAGCCTCAGCTCCCATTCACACCAGAGTCAGATTTACAGAAGCCACCATGTG
ACTGGCGGTTCCTTTGAGGCCCCTCATGCCCCAGGCTCATTTCATGGTGATTATGCCTAT
GGCACATATGCCAGCAATTTCAGCGGTGCCCATGGTTTTCCAGAGTACAGCTACCCTGCA
GACACCTCCTGGCCTGCTGTGGAGCAAGTTCCATCAAGACCGACCTCTCCTGAGAAGTTT
...
...
...
...

I want the output to look like :
>DR-1
1260
>SK-1
1020
>SC-1
1380
>SC-2
960
>TA-1
921
>DR-2
2400
TA-2
>1980
...
...
...
...

I am new to learning python for molecular analysis, so am having trouble with the logic and syntax behind it. I'm looking forward to a warm response from you. Thanks again.

Replace "string protein sequences.fa" with "your_fasta_file.fa" and your file should be at the same path that you created your code.
import Bio.SeqIO as IO
record_dict = IO.to_dict(IO.parse("string_protein_sequences.fa", "fasta"))
for key in record_dict.items():
    print(key[0],"\n ",len(key[1].seq))

Output:
9606.ENSP00000354718
  2817
9606.ENSP00000383118
  191
9606.ENSP00000275364
  381
9606.ENSP00000265563
  404
9606.ENSP00000265717
  418
9606.ENSP00000348461
  211
9606.ENSP00000249071
  192
9606.ENSP00000304283
  192
9606.ENSP00000400175
  193
9606.ENSP00000285735
  193
9606.ENSP00000339467
  191

I hope that's what you need, good luck

Thanks for you quick respone. Unfortunately i didnt have biopython in my windows and also in Linux env. I could not make it run. Can you please give the code to run with general python code?

You can just install it with this command: pip install biopython, if not try the code that is in the second answer

Don't forget to change the file name and write the valid path
with open("string_protein_sequences.fa") as file_content:
    liste=[]
    j=-1
    for seqs in file_content:
        if seqs.startswith('>'):
            j+=1
            print(seqs,lis[j])
            some=0
            i=0
        else:
            i+=1
            some+=len(seqs)
        liste.append(some-i)
    lis=[]
    for j in range(1,len(liste)):
        if liste[j]==0:
            lis.append(liste[j-1])
    lis.append(liste[-1])

Output:
 >9606.ENSP00000354718 (AKAP13)
  2817
 >9606.ENSP00000383118 (CDC42)
  191
 >9606.ENSP00000275364 (GNA12)
  381
 >9606.ENSP00000265563 (PRKAR2A)
  404
 >9606.ENSP00000265717 (PRKAR2B)
  418
 >9606.ENSP00000348461 (RAC1)
  211
 >9606.ENSP00000249071 (RAC2)
  192
 >9606.ENSP00000304283 (RAC3)
  192
 >9606.ENSP00000400175 (RHOA)
  193
 >9606.ENSP00000285735 (RHOC)
  193
 >9606.ENSP00000339467 (RHOG)
  191

you can use a dictionary as well:
fasta = {}
with open("test.fasta", "r") as file:
    for line in file:
        line = line.strip()
        if line[0] == ">":
            header = line
            fasta[header] = ""
        else:
            data = line
            fasta[header] += data

for k in fasta.keys():
    print(f"{k}\n{len(fasta[k])}")

It works great! Thanks. The len information has generated on the bottom row of the header. Can i put this information with the seq id using some name like:

seq1 "len=424"
AATT.........................ATTT

In that case you should use fasta.items() to iterate through the dictionary: for k, v in fasta.items(): print(f'{k} "len={len(fasta[k])}"\n{v}')


text:
Simultaneously get data from multiple applied gates in flowCore
Using the Bioconductor flowCore package, I'm applying two parallel and non-overlapping gates to a gatingSet directly under "root":
library(flowCore)

# Import data
file.name <-system.file("extdata","0877408774.B08",package="flowCore")
x <-read.FCS(file.name, transformation=FALSE)

# Create flowSet
fs <-as(list(x,x), "flowSet")

# Gate
library(flowWorkspace)
gs <- GatingSet(fs)

gate1 <- rectangleGate(filterId = "Gate1", list("FSC-H" = c(300,600), "SSC-H" = c(50,300)))
gs_pop_add(gs, gate1, parent = "root")

gate2 <- rectangleGate(filterId = "Gate2", list("FSC-H" = c(700,800), "SSC-H" = c(500,600)))
gs_pop_add(gs, gate2, parent = "root")

recompute(gs)

gs_get_pop_paths(gs)
# [1] "root" "/Gate1" "/Gate2"

# Plot
library(ggcyto)
autoplot(gs, gs_get_leaf_nodes(gs))

# Extract gated data
gs_pop_get_data(gs, "Gate1") # Works
gs_pop_get_data(gs, "Gate2") # Works
gs_pop_get_data(gs, c("Gate1", "Gate2")) # Does not work
gs_pop_get_data(gs, "Gate1", "Gate2") # Does not work

I want to get the data from the union of gates 1 and 2 but I can't figure out how. The flowCore documentation touches on this topic, but the examples shown are for filters, and I haven't figured out how to get it to work with gates.
If there is a more straightforward way of doing this using another package that would also be of interest, as long as it is compatible with flowSet/flowCore objects.


text:
The biological meaning of the random variables and the responses in Seurat analysis
In the Seurat analysis, if we suppose that Xg and Xr denote the random variables that associate to the expression level of the gene g and of the gene r, respectively. Let Y and Yv represent the responses for Xg and Xr, respectively.
Is it true that biologically, the random variables Xg and Xr represent the distributions of the gene expression measure per cell and the gene expression measure per spot, respectively and the responses Y and Yv represent the range of the random variable Xg and Xr, respectively. The values for Y and Yv are the cell types?
Reference:
https://projecteuclid.org/journals/statistical-science/volume-18/issue-1/Multiple-Hypothesis-Testing-in-Microarray-Experiments/10.1214/ss/1056397487.full
[Method section 2.1]
I try to understand what is the Random variable and the response in this context biologically?
I am confused to make a biological statements about these response. Please help. Thank you.

Its not clear what you mean by "Seurat analysis" Seurat enables the execution of many different analysis types, all of which will be formulated differently in terms of their statistical model and biological interpretation. Further, Seurat performs analyses of single cell RNA seq. The manuscript you have cited is concerned with differential expression testing in bulk microarray data, it is unlikely that the formulation there are helpful.

@IanSudbery I am using the Seurat packages and algorithms and then try to apply the multiple testing analysis referred to the article I cited. So, I got stuck with the Random variables and responses terms they used. I meant Seurat for scRNA-seq datasets.

The article you reference is discussion differential gene expression in bulk samples between conditions. That is, e.g. if take 1 million cells treated with water and 1 million cells treated with a drug and blend each set of 1 million cells together and measure the level of gene expression in each (which will be the average across 1 million cells), how can I tell if gene expression is different after treatment with drug. This is almost certainly not an analysis would want to be doing with scRNA-seq data using Seurat.

@IanSudbery Sorry for misunderstanding, but what I want from the article is to adapt the use of random variables and the term "response" only despite the data they use is bulk samples instead of scRNA-seq datasets. So I am curious the biological meaning of those two  terms.

In the linked article the authors formalize microarray analysis as the study of the joint distributions of $\overrightarrow{X}_i$ and $Y_i$, where $\overrightarrow{X}_i$ is a vector of random variables, the distribution of each of which (i.e. each of $X_{ji}$) is determined by the level of expression of gene j in sample i, and $Y$ is some response or covariate of interest.
The response/covariate is the easiest to interpret in terms of biology. It is the condition. So it might be treated or untreated if we were interested in genes that respond to drug treatment, or it might be diseased vs healthy - this would be covariates. Or it might be "time to death", "body mass index", etc. - these would be responses. The difference between a co-variate and a response is which is the dependent and which is the independent variable. In differential expression we have "gene expression ~  disease" (i.e. X ~ covariate), where ~ is pronounced "is predicted by". In the case of survival analysis "time to death ~ gene expression" (i.e. response ~ X).The malleability of the relationship between X and Y is what is meant by

"The population and sampling mechanism will depend on
the particular application"

A random variable is any variable that we cannot know the value of until it has been measured. Note, this is not "do not" know the value of, its "can not" know the value of.  So for example, in this case of RNA-seq, X is not the gene expression, it is the number of reads. An all-knowing daemon could know the gene expression for a gene in a sample before we measured it, but they couldn't know how many reads we will get because of the randomness introduced by sampling.
Random variables don't represent distributions, they have distributions. That is, if you were to realize many trials of a random variable, the values you obtained would form a distribution.
So the number of reads you get mapping to a particular gene in a sequencing experiment is a random variable. As is the brightness of a spot on a microarray. Its an interesting question to ask if $Y$ is a random variable. The answer is "it depends" - and what it depends on is how good we are at controlling an expriment. You might argue that "drug treatment" is not a random variable because we treat the cells or not - we know, but that, say "time to death" or even "age" is a random variable, because we don't determine it ahead of the experiment. Either way, this subtlety is covered by the above quoted sentence.
So, in a traditional bulk sequencing experiment we have a matrix of random variables associated with gene expression $X_{ji}$ - the number of reads for gene j in sample i, and $Y_i$ the condition for sample i. Note that $Y$ is indexed by the sample, not the gene - every gene in the sample uses the same value of $Y$. And for each gene j we test the model X_i ~ Y_i - the expression in sample i is predicted by the condition of sample i.
In a single cell sequencing our matrix of random variables associated with gene expression is actaully a tensor, that is, it has three indexes, not two: $X_{jki}$ - that is the number of reads from gene j in cell k in sample i. In a scarily large number of single cell sequencing experiments, i only has one value (i.e. there is only one sample). And in most single cell experiments I've seen $Y_i$ only takes one value, even if there are multiple samples. Even if $Y$ takes multiple values, its not clear what sort of model you'd fit to the structure X_{ki} ~ Y_i - i.e predicting a matrix of many cells by a single vector of covariate values.
I think one analysis you can do with Seurat is to treat each cell as a separate sample, and each cluster as a condition, and thus try to find gene DE between cell types. I'm not convinced how great this sort of analysis, but here, the i in $X_{ji}$ would represent which cell a measurement was from and $Y_{i}$ would represent which cluster. But this sort of analysis seems prone to all sorts of date-peeking and training set leakage problems.

if we suppose that Xg and Xr denote the random variables that associate to the expression level of the gene g and of the gene r, respectively. Let Y and Yv represent the responses for Xg and Xr, respectively.

First, $Y$ is a vector of which $Y_v$ is a component, so $Y$ contains $Y_v$. And in the formulation presented in the linked paper (and in almost all RNAseq analysis literature), the response is the same for two different genes, so there is only one $Y$ across both genes g and r.

the random variables Xg and Xr represent the distributions of the gene expression measure per cell and the gene expression measure per spot, respectively

I'm not quite sure what you mean by "per spot". Spots are a microarray term, and don't really have any counterpart in RNA-seq analysis. They were, literally, the spots of DNA oligos printed onto the glass slides used in microarray analysis.
$x_{gi}$ is the number of reads mapped to gene g in cell i, and has a distribution determined by the level of expression of gene g in the cell i.
$x_{ri}$ is the number of reads mapped to gene r in cell i, and has a distribution determined by the level of expression of gene r in cell i.

and the responses Y and Yv represent the range of the random variable Xg and Xr, respectively

No. The response $Y$ represents some biological dicotomy you would like to associate with gene expression.  You could use cell type, but this is neither response, nor covariate, it is not measured independently of X, but is simply a phantom of the correlation structure in $X_{ji}$, so I can see why the "response" type language would be confusing.

Thank you very much. The terms "response" and "covariate" are now clear to me. Also, your explanations are really useful. I do appreciate it. I accept this answer


text:
Hg38 annotation tracks retrieval
I want to retrieve an annotation object (as GRanges) containing different UCSC annotation tracks such as CpG Islands, TSS, TFBS and promoters for hg38 Human reference genome. But "UCSC Genes" was phased out for hg38. What could be the best resource for this ?
Thanks
RT

The UCSC Data Integrator can combine multiple tracks (CpG Islands, TFBS, cCREs, etc.) into one output annotation file, and it works with any assembly you choose. You can select which output fields you want, then download the tab-delimited file.
From there, you can use R, Python, and/or GNU tools such as awk and sed to coerce the data into a format which is useful for your purposes. You can then load it into R as a data frame and use the columns as arguments in the GRanges constructor to create a GRanges object.


text:
How to calculate RMSD and binding energy from the already docked ligand-receptor
I have the following PDB file that looks like this.
It is downloadable here. This PDB already contains ligand and receptor together.
The red molecule is the ligand and the yellow is the receptor.
They're already docked with ClusPro.
My question is how to calculate the RMSD and binding energy from it?
Any existing tools (with or without command line) that can help?
I don't have any data to try it out at hand but have a look at the documentation I linked. There you can find ways of constructing an Ensemble which can then be superposed and RMSDs calculated, even pairwise.
Reference is the red one (receptor). I have hundreds of other complexed ligands. I want to calculate how good are those complexes. What's the reasonable measure to do that?
Example code using Prody? please :-)

Could you figure it out how to calculate the binding energy? which software did you use? Thanks,

NB. This is a long comment as opposed to an answer which requires a colossal working chunk of code.

Control
Science is all about controls failing. A great control for docking is the barnase-barnstar complex (PDB:1BRS): both small, but form a strong interactions. If this is done, the RMSD against the original structure is very meaningful and useful.
RMSF
A RMSD is a metric of how much two structures differ. In an ensemble, the metric of how much a structure and the ensemble average is called a RMSF (root-mean square fluctuation). Here the RMSF is sought as the RMSD from first/top pose or worse the manually thrown together pre-docking pose does not mean much.
Clusters
The decoys (PBDs trails) need to be clustered into different ensembles first, before a RMSF can be calculated.
∆∆G and water
For the calculation of difference in Gibbs free energy between the docked minus the docked poses, the structures need to be energy minimised with the forcefield scorefunction that is to be used for scoring.
Water is very important for protein-protein interactions, but adds a lot of complication. Say Gromacs is used, gmx energy will return the single state (conformer) Gibbs free energy, but it requires a water box. For the difference, the water box of undocked needs to be the same size and preferably with the same amount of molecules. A much simpler solution is using an implicit solvent model, e.g. Rosetta. As a bonus, the energy minimisation (relax) is relative fast. And specific explicit waters can be added. However, calculating ∆∆G for every decoy will take hundreds of CPU-hours and is a bit nonsensical.
A better approach is to look at residue conservation (consurf), cross-linking mass spec data if available, and if human post-translational modifications (phosphosite) and healthy variants (gnomAD). But ultimately wet lab validation by ITC of different mutants/labelling based upon the docking results etc. is required.


text:
Is New Tuxedo approach really better than the earlier ones?
I am planning to perform multiple routine NGS analyses for a hospital lab on patients' samples on an everyday basis.
Before I used a new tuxedo approach with ballgown etc and was rather satisfied with it.
I wonder if for the above-mentioned job it's ok to use a ballgown or, perhaps, that is not the best choice?

I'm fairly used to analyse NGS illumina data, but what is the New Tuxedo approach you mention? How do you measure best choice?

Which kind of data are you working with? Doesn't a clinical setting require some kind of legal approval and validation of the pipelines you use? Please specify which data you have and what even the analysis goal is.

Thanks I am talking about this approach https://www.nature.com/articles/nprot.2016.095


text:
Pseudotime analysis on multiple samples
I am using Monocle3 to run pseudo time analysis on samples at two different conditions (samples are paired). My question is, is it valid to combine all the cells from samples at specific condition when clustering, learning graph and ordering cells? or should I pick up only one sample per condition and run the analysis on?
My understanding is that pseudotime analysis uses the expression data of cells to try to infer the trajectory of cells, so I believe it might not be correct when I pool all the cells that might not have a true trajectory?
Thanks

I would try it both ways (i.e. run pseudotime analysis on individual and combined samples). If you have the same transcriptional trajectories in multiple samples then you might be able to capture this via combining samples, and arguably the most interesting signatures are generalizable across samples. On the other hand, there might be batch effects that disrupt the signal of a combined analysis.
I think I understand your concern is that different samples have independent cells and so they should not form part of one trajectory. However, even within a sample, we are not capturing the true trajectory as this would require sampling the same individual cell multiple times at different time points, which is not currently possible because scRNA-seq requires lysing the cells (see also live-seq).

That makes sense. I found a similar framework with more formal analysis if we should or shouldn't use a common trajectory (but between different conditions here) in this paper (https://www.biorxiv.org/content/10.1101/2021.03.09.433671v1). Anyway, thanks!


text:
Clustering by individual tissue population post integration of two datasets
Sorry if this question has been asked. Post batch correction and integration of two datasets, where cells of multiple different tissue types are present, I want to run clustering on cells of corresponding tissue types between the two datasets. Should I simply subset the integrated object to get each tissue type, scale the data and run PCA and then run clustering on it?
Thanks so much for the clarifications.


text:
use Kallisto in galaxy
I want to use Kallisto for sequence alignment in Galaxy. Its description is:

a program for quantifying abundances of transcripts from bulk and
single-cell RNA-Seq data, or more generally of target sequences using
high-throughput sequencing reads.

However, I found this field empty:

there is no option available for the reference transcriptome.

How I use a reference sequence?
Can I leave the other configuration by default?

this is the link of the tool:  https://usegalaxy.org/root?tool_id=toolshed.g2.bx.psu.edu/repos/iuc/kallisto_quant/kallisto_quant/0.43.1.3

You have to click on the dropdown menu and select "Use transcriptome from history", which will allow you to specify your own file.

I suspect the default will be built-in human transcriptome.
Solution is not tested, so let us know if it actually works.


text:
A database for RefSeq protein accession IDs
Are there any databases or tools that contain a list of all RefSeq protein IDs for the interest organism/microorganism, e.g. get a specific microorganism name and return all RefSeq Protein IDs that have been identified?

Please consider accepting the answer below - its good for site stats and its a good answer IMO

How about using RefSeq itself? Say, you are interested in all E. coli proteins, you can use the query "Escherichia coli"[Organism] AND RefSeq[Filter] here.
If you know the Taxonomy identifier for the specific organism of interest to you, it can also be used. For example, if you are interested in only "Escherichia coli str. K-12 substr. DH10B", you can use the taxonomy identifier 316385 (obtained from here) in the query txid316385[Organism] AND RefSeq[Filter].
To extract a list of accessions and save them to a file on disk, use the "Send To" menu at the top right corner, choose "File" as destination followed by "Accession List" as the format and download the file. Alternatively, you can use E-utilities API or the EntrezDirect tool on the Unix command line to download these data programmatically.
To address the question about E. coli proteins shared between different species... There's no need to parse FTP files to get all the RefSeq proteins for a bacteria species. While it's true that some RefSeq WP proteins have taxids above species level, like Enterobacterales WP_003861789.1 where the identical protein is found in genomes from multiple species of the order Enterobacterales, those proteins are still retrieved by a species-level organism search. NCBI implemented some extra handling for the RefSeq WP proteins to treat the MULTISPECIES WPs like they belong to all of the species with genomes on which they're annotated. So of the 56947 WPs assigned as "Bacteria", 7401 are found by "Escherichia coli"[Organism], and the species-level search will give you the full set of proteins.
You could also use the AutonomousProtein2Genomic.gz file from the RefSeq FTP path which has a mapping of how all WPs map to specific genomic sequences, but it's a really big file and it would take a long time to grep for "Escherichia coli".

Thanks warmly for the response.
Actually, I was looking for a database in which I can extract all RefSeq protein IDs as a .csv or .xls file.

I would say,suggested query is not working as expected for refseq.Some escherichia coli refseq proteins have “lower” taxonomy rank ,for example “”Escherichia- proteins which shared between different Escherichia species.Technically ,I know only one way for now extract ALL Escherichia Coli refseq proteins

retrieve all corresponding refseq assemblies from ftp and parse them


text:
"variable.features.n" in SCTransform
Whats the role of "variable.features.n" in SCTransform function of Seurat? I try setting "variable.features.n" at three different numbers (2000, 5000, 10000). In the downstream cells are clustering in different groups?

variable.features.n sets the number of features (you can think of genes in the case of scRNA-seq) you would like to use for the downstream steps such as clustering. Basically, using a gene that is expressed more or less at similar levels across different cell types would not be informative in terms of differentiating (for example via clustering) these cell types, this parameter helps you to choose the "most informative" (most variable across the whole data set) features.


text:
Making colnames by number and character
I have a very big matrix I want to put colnames as cell1, cell2, et,.. but I get error
> colnames(new.dat)
NULL
> View(new.dat)
> dim(new.dat)
[1]     34 263086
> colnames(new.dat)=make.names(c(1: 263086, "cell"))
Error in dimnames(x) <- dn :
  length of 'dimnames' [2] not equal to array extent
> colnames(new.dat) <- paste0("cell", seq_along(new.dat))
Error in dimnames(x) <- dn :
  length of 'dimnames' [2] not equal to array extent
>

any help?

You should really test your code:
make.names(c(1: 5, "cell"))
[1] "X1"   "X2"   "X3"   "X4"   "X5"   "cell"

This is obviously wrong.
Better:
colnames(new.dat) <- paste0("cell", 1:ncol(new.dat))


text:
Pulling a numbered chromosome range file given a gene location from a lookup table ideally from command line or R
have a folder with roughly 1000 vcf files which have divided the human genome into chunks, the folder looks like this:
- main_programme_aggregated_chr1_1_194789237.vcf.gz
....
- main_programme_aggregated_chrX_14235_8759493845.vcf.gz
etc

I will divide this file name into:
*_chr[1-22, X,Y]_A_B.*

I have a separate tab divided file with the names of genes and their genomic coordinates (all the genes in the human genome on build 38) - file1:
chr1 19853 90835 WASH7P ENSG0000022454
...
chr10 38732 390853 JHBF ENSG0000578382
..
chrX 9532 908032 LOTY ENSG00005847

etc
I would like to use the gene name from the 4th column of file1 to find the relevant chunked file it falls within. So I would need to take column 1 of file 1, match it to chr[1-22, X or Y] and then find the chunk were the range of columns2+3 fall between (A-B) in the file name.
I do not need to go into the files themselves, just use their names.
I have been just manually looking up the gene coordinates and the relevant chunk but am sure this can be automated.
An example would be that I want to get the chunk file that the gene WASH7P falls under, so the relevant file would be main_programme_aggregated_chr1_1_194789237.vcf.gz.
Ideally I could do this from the command line and the relevant chunk file then be inputted as a variable into some downstream workflows I have.
Would anyone be able to help? Many thanks for your time

Here is a solution with bash:
#!/bin/bash

# 1 - create an array of coordinate files

# use nullglob in case there are no matching files,
# not that important for this toy case
shopt -s nullglob

# make a(n) array/list of vcf file in the current (./) directory
coordinate_files=(./*.vcf.gz)

# going through each line of "file1"
while read line; do

    # 2 - get the gene related info from "file1"
    # 2.1 - change line to an array (line split by whitespace)
    line_array=( $line )
    # 2.2 - get the 4th column, gene name and chr, the first column
    gene=${line_array[3]}
    chr=${line_array[0]}
    # 2.3 get coordinates
    start=${line_array[1]}
    end=${line_array[2]}

    # loop through each file in the array created at step 1
    for ((i=0; i<${#coordinate_files[@]}; i++)); do

        # parse the info hardcoded in the filename
        chr_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $4}')
        start_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $5}')
        end_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $6}')
        end_in_filename=${end_in_filename%.vcf.gz}

        if [ "$chr" == "$chr_in_filename" ]
        then
            if [ "$start" -ge "$start_in_filename" ] && [ "$end" -le "$end_in_filename" ]
            then
                echo $gene
                echo ${coordinate_files[$i]}
            fi
        fi
    done

done <file1

The input is the two example coordinate files + "file1" in the question and the output is:
WASH7P
chr1 19853 90835 WASH7P ENSG0000022454
./main_programme_aggregated_chr1_1_194789237.vcf.gz

Many thanks for taking the time to answer. I understand the principle behind this but the actual file chunks look like this: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz (as an example), how would I adjust your script for this?

The script looks for underscores in the filenames of the coordinate files, so $4 becomes the chromosome,chr_in_filename: so for main_programme_aggregated_chr1_1_194789237.vcf.gz, the chromosome is chr1 . You can play with the numbers in the for loop. But since the number of "fields" in the filenames in your original post and what you pasted in your comment is the same, the code should work fine. How do you run this code and what is the output?

Below is a Python solution:
import os

files = {}

for root, dirs, files in os.walk('/path/to/vcf/files'):
    for file in files:
        # Example: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz
        chrom = file.split('_')[3]
        start = int(file.split('_')[4])
        end = int(file.split('_')[5])
        files[os.path.abspath(file)] = {'chrom': chrom, 'start': start, 'end': end}
    break

genes = {}

with open(file1) as f:
    for line in f:
        fields = line.strip().split('\t')
        chrom = fields[0]
        start = int(fields[1])
        end = int(fields[2])
        gene = fields[3]
        genes[gene] = {'chrom': chrom, 'start': start, 'end': end}

mapping = {}

for gene in genes:
    for file in files:
        if (genes[gene]['chrom'] == files[file]['chrom'] and
            genes[gene]['start'] >= files[file]['start'] and
            genes[gene]['end'] <= files[file]['end']):
            mapping[gene] = file
            break

with open('mappings.txt', 'w') as f:
    for gene in mapping:
        f.write(gene + '\t' + file + '\n')

Many thanks for taking the time to answer. I understand the principle behind this but the actual file chunks look like this: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz (as an example), how would I adjust your script for this?

@tacrolimus, if all the files look similar to your example (main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz), then above script should run just fine. Basically, the script will split each file's name by the underscore character ('_') and then picks the 3rd item as chromosome, 4th item as start, and 5th item as end (0-based indexing). Let me know in the comment if somehow it fails.

Worked well thanks! I have accepted the shell one as that my more "native" language but really appreciate you taking the time to answer and has spurned me to look into python more!


text:
How do I generate sequence logo of ATGC against each position in R?
I am new in R sorry for the naive question. I want to generate the seqlogo for each position using R. Can anyone guide me? Thank you so much in advance.
1  GATGGGGTGGGCAACTTAGAACTT
2  ATAACTGT
3  CAGAGGCGAACT
4  GGGCTAGGG
5  GGGGGGTGGCAGCACTAGC
6  GGGTAATATG
7  GTCGAGCGGT

Best

It's good practice to google before asking on online forums. Searching for "r seqlogo" pointed me to the seqLogo BioConductor package, which might be really helpful to you: https://bioconductor.org/packages/release/bioc/vignettes/seqLogo/inst/doc/seqLogo.html
You will need your sequences aligned though.


text:
Help with DIA-Mass Spectrometry data analysis with several conditions (limma?)
I am trying to learn how to analyse normalised DIA-MS data and I am struggling with it ://
The original dataset I got is (6 conditions (2 samples each)) with 3 technical replicates (total: 36 sample columns x 1128 proteins in rows). I removed the proteins with more than 80% of NA, and impute with kNN the rest of missing values.
Then I collapse the replicates doing the average of the 3 replicates (is this correct?)
So now, I have a dataframe with 12 columns x 1072 rows.
The question is... I am not sure how to perform the differential analysis from here (assuming I did well the previous steps lol)
Could someone help me with the tools or code to compare all the 6 conditions between them?
I read something about Limma package, but not sure if it can be used in this DIA-ms data and neither sure of how to use limma for it (I already check the vignette with no succeed)
UPDATE: I tried to use Limma as follows:
fit1 <- lmFit(norm.data_collapsed, design)
fit2 <- contrasts.fit(fit1, contrasts = contrast)
fit3 <- eBayes(fit2)
limma.res.01 <- topTable(fit3, coef = comb.contrast[1], n= Inf)

This returns me a non-sense values for logFC or significance:

Any help would be very welcome,
Best regards

Regarding the big FC values, I just realise the Limma lmfit function requires as input a log2 matrix... so I just transformed my data and now makes more sense. Still I have the doubt regarding the collapsing technical replicates and and if some one can verify the pipeline would be like that.
Thanks!!!
Fer

Please could you mark this answer as accepted? That would be very welcome.


text:
What are some ways to error correct Oxford Nanopore long read sequencing?
I am sequencing long read genomic sequences to assemble MHC region haplotypes in a non model organism using Cas9 Sequencing Kit (SQK-CS9109) using flongle adaptors in the minION from Oxford Nanopore Technologies (ONT). However the error rate in genomic sequences results from ONT has been over 14%.
I know that there is error correction using Illumina sequencing and I was hoping anyone knows more about how this is accomplished?
Additionally, there are bioinformatic tools to error correct these errors? Has anyone accomplished diminishing this big error rate using bioinformatic tools?

How are you evaluating sequencing error rate? My most recent re-calls of 2017 sequences are demonstrating median single-read accuracy over 96%.
Before considering Illumina, it'd be worth it to do an initial correction using nanopore-only reads. This will make sure that the best results are obtained from the nanopore signal. First, make sure that the reads have been recalled using the latest highest-accuracy basecaller (currently guppy v5.0.11):
guppy_basecaller -c dna_r9.4.1_450bps_sup.cfg ... # for 9.4.1 reads
guppy_basecaller -c dna_r10.3_450bps_sup.cfg ... # for 10.3 reads

If you don't have any reference sequences available, reads can be corrected using the other reads by using canu:
canu -correct -p nonModel -d nonModel genomeSize=10m -nanopore reads.fastq

If you do have reference sequences available, then reads can be corrected using megalodon. I don't have any experience with megalodon, but have seen a couple of papers demonstrating that it has good accuracy, especially for methylation calling.
After all that, if you're still not happy with the accuracy, you can use Pilon for read correction using the Illumina reads. Pilon is used for correcting a reference sequence, so you'll need to create your haplotype assembly first before using Pilon. Even when using Pilon, I'd still recommend only using the Illumina reads for local INDEL correction. Pileup nanopore reads are very accurate for substitutions, especially when considering [lack of] strand bias when evaluating variants.
Update:
Given that you're interested in determining haplotypes across the MHC, which is something that performs poorly for short reads even for well-annotated genomes, I think you're better of sticking with reference-free correction, and categorising haplotypes while taking taking this into account. If you just need to be able to distinguish haplotypes, I recommend doing homopolymer compression on the resultant sequences, as hompolymer inaccuracies are the most common error mode for nanopore reads. Even if that doesn't get you all the way to perfect haplotypes, it might get close enough that something like CDHIT can be used for binning similar sequences.
HLA typing in humans can be done from nanopore reads (e.g. see here and here), but don't expect perfect base-level accuracy. What I'm suggesting is something similar to a typing process, but without an established reference, i.e. binning reads based on their similarity to each other.

Thank you very much for this amazing response it is very comprehensive. What do you think are the major issues with haplotype assembly when sequencing the MHC as we are interested in this region? While there is a reference genome for this non model organism it is poorly annotated (has not been updated for 10 years) and the fact that we are interested in the MHC region in particular makes things even more complex. Would you have any additional recommendations about error correction in the MHC region due to the high polymorphism?

@mk894 I haven't seen papers or preprints that can get good haplotype assembly in HLA only using nanopore reads from one sample. A major problem is that nanopore error correctors can't reliably distinguish sequencing errors from real polymorphisms. This may lead to false correction and mess up downstream assembly.


text:
Intrinsically disordered protein
This question has also been asked on Biostars
The position that I have to visualize is in the 2700 position and this position exists in two ranges in https://www.uniprot.org/uniprot/Q12802#family_and_domains the first range is Interaction with ESR1 and the second region is disordered, the first question is which range I should be based on? and the second is Why intrinsically disordered protein(IDP) can't be visualized in 3D?

Asking how to visualise the disordered region of a protein is a bit like asking how to visualise the location of an electron.
An intrinsically disordered protein, by definition, has regions that are disordered. In other words, their location is variable under observation, so they cannot be precisely placed in a 3D model.

Thank you for the answer @gringer, I understood.


text:
Comparing the genetic similarity of human genome and mice/monkey genome - where can I find a list of genes?
I am trying to start a side project where I verify the claim that human and mice/monkey share 99% of their genes. I hope to learn the basics of doing genomic analysis.
I would like guidance about where I can find the sequences I need, and programming languages or applications that I should use. For the sequences, I would like for example, a list of all genes for each species. In the case that the programming language or application is foreign, some basic guidance would be helpful.
I have programming experience with matlab and python, and have learnt R previously, though forgotten. Would be able to learn a new language.
Alternatively, I think this may be a very large subject. If there is a link to a tutorial/lectures it will also be helpful.
Thanks in advance.

Thanks for including context around your question. Remember to keep your posts to a single question. If you have a series of related or linked questions, please create one post for each question.

As you say, this is a big project, but I like the idea of having a themed set of questions around a single workflow; I think it's a great way to attract interest and keep people coming back to this site.

For downloading lists of genes, together with associated features, I like using Ensembl Biomart:
http://ensembl.org/biomart/martview/
In this case, you can "CHOOSE DATABASE -> Ensembl Genes", then "CHOOSE DATASET -> Human Genes" or "CHOOSE DATASET -> Chimpanzee genes" to get to a table selection.
Clicking on "Attributes" on the left hand sidebar brings you to a selector for different fields to add to the table.
Clicking on the "Results" button at the top left brings you to a form for choosing how results should be output.
Bear in mind that gene names are not likely to match. Finding the intersection of gene sets between species is a complicated process, usually requiring knowledge about gene copies, paralogs, orthologs, and sequence comparison.
It's also likely that the human genes are more annotated than the chimpanzee genes, and there might be genes present in the human list that are also in the chimpanzee genome, but have not yet been discovered.

Thanks, just what I need


text:
How do I get GO annotations for a list of UniProt IDs?
I have a list of UniProt ids that I want to get Gene Ontology annotations for. I need this information because I want this high-level information as an input to a neural network. The model I wish to develop is inspired by this paper: https://www.ijcai.org/Proceedings/2018/0468.pdf
Can anyone advise how I can programmatically access this information in the
Gene Ontology Annotation (GOA) Database? https://www.ebi.ac.uk/GOA/

The link to the FTP for the GOA database files is listed on the GOA Downloads page. The file containing the mapping info you seek, goa_uniprot.all, comes in two formats, .gaf and .gpa. The README in the FTP directory linked provides details on the structure of the files.
Below is the first 20 lines of the .gaf file after uncompressing. Column 2 contains the UniProtID and column 5 contains the associated GO term. (Note that column 4 is empty in the snippet below)
>$ head goa_uniprot_all.gaf  -n 20
!gaf-version: 2.2
!
!The set of protein accessions included in this file is based on UniProt reference proteomes, which provide one protein per gene.
!They include the protein sequences annotated in Swiss-Prot or the longest TrEMBL transcript if there is no Swiss-Prot record.
!
!date-generated: 2021-02-12 21:14
!generated-by: UniProt
!go-version: http://purl.obolibrary.org/obo/go/releases/2021-02-08/extensions/go-plus.owl
!
UniProtKB   A0A009GUA7  J508_4179       GO:0003677  GO_REF:0000002  IEA InterPro:IPR006119|InterPro:IPR036162   F   Resolvase/invertase-type recombinase catalytic domain-containing protein    J508_4179   protein taxon:1310609   20201128    InterPro
UniProtKB   A0A009GUA7  J508_4179       GO:0006310  GO_REF:0000002  IEA InterPro:IPR006119|InterPro:IPR036162   P   Resolvase/invertase-type recombinase catalytic domain-containing protein    J508_4179   protein taxon:1310609   20201128    InterPro
UniProtKB   A0A009GUA7  J508_4179       GO:0000150  GO_REF:0000002  IEA InterPro:IPR006119|InterPro:IPR036162   F   Resolvase/invertase-type recombinase catalytic domain-containing protein    J508_4179   protein taxon:1310609   20201128    InterPro
UniProtKB   A0A009HCR2  J517_0313       GO:0005886  GO_REF:0000044  IEA UniProtKB-SubCell:SL-0039   C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0016020  GO_REF:0000043  IEA UniProtKB-KW:KW-0472    C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0016021  GO_REF:0000043  IEA UniProtKB-KW:KW-0812    C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0015031  GO_REF:0000043  IEA UniProtKB-KW:KW-0653    P   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0005886  GO_REF:0000043  IEA UniProtKB-KW:KW-1003    C   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    UniProt
UniProtKB   A0A009HCR2  J517_0313       GO:0055085  GO_REF:0000002  IEA InterPro:IPR003400  P   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    InterPro
UniProtKB   A0A009HCR2  J517_0313       GO:0022857  GO_REF:0000002  IEA InterPro:IPR003400  F   Biopolymer transport ExbD/TolR family protein   J517_0313   protein taxon:1310618   20201128    InterPro
UniProtKB   A0A009HLK9  J512_2607       GO:0006355  GO_REF:0000002  IEA InterPro:IPR000847|InterPro:IPR037402   P   HTH lysR-type domain-containing protein J512_2607   protein taxon:1310613   20201128    InterPro

That is exactly the information that I need. I am using a Mac and I have used The Unarchiver and Archive Utility to try and unzip the file. I get a message saying that the file is incomplete. How did you unzip the goa_uniprot_all.gaf file?

@ChemBot that isn't a zip file, you can't unzip it. It is a gzipped file, so you need gzip. Open a terminal and run gunzip goa_uniprot.all.gz.


text:
Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'
I am trying to open .bam files in C++ to extract reads occurring at specific scaffolds and loci. I essentially want to call "samtools view sample.bam -o sample.sam scaffold:pos-pos" from C++. I have tried system("samtools view sample.bam -o sample.sam scaffold:pos-pos") and can explain the issues if you want to know, but for the sake of conciseness I will leave the details out of this post. So, I am using SeqAn. When I run my script, I get the following error:
terminate called after throwing an instance of 'seqan::UnknownExtensionError'
  what():  Unknown file extension of copyPolly_SRR6511930.bam: iostream error

I am running my program on my institution's supercomputer clusters (Texas A&M Terra cluster), to which I do not have sudo access. I have loaded the following modules:
GCC/9.3.0
SeqAn/2.4.0-GCCcore-9.3.0
zlib/1.2.11-GCCcore-9.3.0

Here is my minimal script:
#include <seqan/bam_io.h>
#include <zlib.h>

using namespace seqan;

int main(){
    CharString bamFileName = "copyPolly_SRR6511930.bam";
// Open input file.
    BamFileIn bamFileIn(toCString(bamFileName));
    if (!open(bamFileIn, toCString(bamFileName))){
        std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
        return 1;
    }

    unsigned numUnmappedReads = 0;

    try{
        // Read header.
        BamHeader header;
        readHeader(header, bamFileIn);

        // Read records.
        BamAlignmentRecord record;
        while (!atEnd(bamFileIn)){
            readRecord(record, bamFileIn);
            if (hasFlagUnmapped(record))
                numUnmappedReads += 1;
        }
    }
    catch (Exception const & e)
    {
        std::cout << "ERROR: " << e.what() << std::endl;
        return 1;
    }
    return 0;
}

I have tried the following solutions:
-Adding #define SEQAN_HAS_ZLIB (1) to my script's header
-Adding -D SEQAN_HAS_ZLIB to my compilation such that it reads g++ script.cpp -D SEQAN_HAS_ZLIB
-Adding -lz to my compilation such that it reads g++ script.cpp -lz
I am not dead-set on using SeqAn for this project, so if anyone knows of another C++ library which  extracts specific reads from .bam files, please share.
Update: I've never used HTSlib and it seems to have a fairly steep learning curve, so if anyone has an answer for the SeqAn error, that would still be very much appreciated.

The helpdesk for my institution's supercomputers showed me the following flags:
g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread

This fixed the error. I'm now running into different errors and am suspicious that I am getting them because I am using SeqAn2 rather than SeqAn3. I'm going to try HTSlib.

Htslib is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM:
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  while (sam_read1(fp, h, b) >= 0) {
    if (b->core.tid < 0) continue;
    int end_pos = bam_endpos(b);
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos, end_pos);
  }
  bam_destroy1(b);
  bam_hdr_destroy(h);
  sam_close(fp);
}

PS: the above started as a comment for "relevant but minor or transient information to" the original question. I didn't consider htslib a direct answer to the question, but because this had been changed to an answer, I added some details anyway. By the way, seqan2 has largely been discontinued. Use seqan3.

Thank you for this! I tested this script and it worked. I will try to edit it to suit my needs. Also, I am using SeqAn2 because SeqAn3 is not available on my institution's supercomputers. If HTSlib doesn't pan out or if no one can figure this error out, I will email them to request that it be installed.

I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the SeqAn BAM tutorial, it looks like you need to use the BamFileIn object for reading BAM/SAM files:
#include <seqan/bam_io.h>
using namespace seqan;

int main(){
    CharString bamFileName = "copyPolly_SRR6511930.bam";
    // Open input file, BamFileIn can read SAM and BAM files.
    BamFileIn bamFileIn;
    if (!open(bamFileIn, toCString(bamFileName))) {
        std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
        return 1;
    }
    return 0;
}

Edit: after seeing the updated code... I've got nothing, sorry. It all seems like reasonable code that shouldn't be generating errors.

Sorry about that, I did not include all of my code. I have updated my post with the rest of the code which led to my error.

I think SEQAN_HAS_ZLIB flag not works correctly as it supposed to be.
So you can check the seqan library to see if SEQAN_HAS_ZLIB works correctly, and seek the reason of extension error.
Seqan library is a mature library which using template to define general-purpose classes and structs.
1.The reason for extension error.
Step into the seqan library:
When you call BamfileIn: BamFileIn bamFileIn(toCString(bamFileName));
An object of BamFileIn struct is created, and bamFileName file is opened by the struct construction function.
And when this object is created, some extra information is wrapped with the object to check the type of file.
This information include suffix categories. In bam_file.h header file.
template <typename TDirection, typename TSpec>
struct FileFormat<FormattedFile<Bam, TDirection, TSpec> >
{
#if SEQAN_HAS_ZLIB
    typedef TagSelector<
                TagList<Bam,
                TagList<Sam
                > >
            > Type;
#else
    typedef Sam Type;
#endif
};

As you can see, if SEQAN_HAS_ZLIB is setted, "Bam" file will be added to TagSelector which will be wrapped into the BamFileIn struct and redefined by typedef as format.
Then _open() function is called to open Bam file, and extension is checked in guess_format.h
template <typename TFilename, typename TTagList>
inline bool guessFormatFromFilename(TFilename const &fname, TagSelector<TTagList> &format)
{
    typedef typename TTagList::Type TFormat;

    if (value(format) == -1 || isEqual(format, TFormat()))
    {
        // if tagId is set to -1 (auto-detect) or the current format (TFormatTag) then test for TFormatTag format
        if (guessFormatFromFilename(fname, TFormat()))
        {
            assign(format, TFormat());
            return true;
        }
    }
    return guessFormatFromFilename(fname, static_cast<typename TagSelector<TTagList>::Base &>(format));
}

template <typename TFilename, typename TFormat_>
inline bool guessFormatFromFilename(TFilename const & fileName, Tag<TFormat_> /*formatTag*/)
{
    typedef typename Value<TFilename>::Type                                     TValue;
    typedef ModifiedString<TFilename const, ModView<FunctorLowcase<TValue> > >    TLowcase;
    typedef Tag<TFormat_>                                                       TFormat;

    TLowcase lowcaseFileName(fileName);
    for (unsigned i = 0; i < sizeof(FileExtensions<TFormat>::VALUE) / sizeof(char*); ++i)
        if (endsWith(lowcaseFileName, lowerString(FileExtensions<TFormat>::VALUE[i])))
            return true;

    return false;
}

The above code is two instance of  guessFormatFromFilename function.
It's pretty complex here, but if SEQAN_HAS_ZLIB is set correctly, TagSelector will include bam and sam, and the first guessFormatFromFilename above will return true in if condition, and the error will not happend. But if SEQAN_HAS_ZLIB is not working, second formt guessFormatFromFilename function will be called and only "sam" suffix is exist in TFormat, and it will be exhibited by FileExtensions<TFormat>::VALUE[i]. You can using std::cerr to print the suffix on your screen to see if there is only "sam" suffix is exist in TFormat.
Althrough you add -DSEAN_HAS_ZLIB, I think you could check the code in library to see if it really works fine.
If SEQAN_HAS_ZLIB is working, you could check other variables or functiones related with open bam process.These function are in bam_file.h, formatted_file.h, fundamental_tag.h and so on.


text:
Is there a way to calculate the Ka/Ks value in R/unix?
I am investigating the rate of synonymous non synonymous substitutions and would like to find the Ka/Ks value also known as dn and ds between multiple alleles? I want to investigate multiple alleles because I am examining the MHC region which is highly polymorphic. Any suggestions would be very helpful!

I have also used the seqinr package for Ka/Ks

A while back, I used PAML for this purpose for a phylogenetics project. It has Windows and macOS executables but for *nix, you will need to compile it from source code.


text:
Does it make sense to perform pathway enrichment analysis on a list of mutated genes from a VCF?
I have a list of mutated genes from a VCF. Does it make sense to perform a pathway enrichment analysis using PANTHER or STRING as examples to know which pathway is affected?
EDIT: to reply to @terdon
We sequence mostly human genome data using Illumina sequencers. As for the type of data, we perform mostly Whole-Exome Sequencing (WES), Whole-Genome Sequencing (WGS).
At the end, after performing all pre-processing, we obtain the VCF files from HaplotypeCaller using GATK. Then we proceed to annotate using VEP.

That will depend entirely on what your data are. Did you sequence an entire genome? An entire exome? A small panel of target genes? Have you separated pathogenic from non pathogenic variants? How about compound heterozygotes? Please [edit] your question and give us some context to understand what you are trying to do. "A list of mutated genes" doesn't really tell us anything.

@terdon How and why does it matter what type of data is it coming from? Whether it's from an entire genome or exome panel? I can understand the question separating the pathogenic from non-pathogenic, but a mutation is a mutation whether coming into exons or introns, right? By compound heterozygotes, do you mean a multi-allelic mutation of a variant?

It matters because you want to do an enrichment analysis. If you're just looking at a small panel of genes, that would be a completely different proposition to looking at WGS or even WES data. If you are working with WGS data and have a set of affected samples and a set of controls, then enrichment analysis could make sense (but that would again depend on what you are trying to test for). However, you should first identify the variants (or genes with variants) that segregate with your affected samples. And by compound het I mean cases where you have two het variants in the same gene [cont...]

... and it is the combination of the two variants that results in the affected phenotype, while neither of the variants would be sufficient by itself.

Thank you. There's still something that I didn't understand completely. You mentioned that if I have WGS data, it would make sense to perform enrichment analysis of affected vs controls. Why would I need a set of control when I want to know simply which pathway is mutated by enriching the mutated genes by comparing them to all genes? For example, the patient I am analyzing is diabetic. I expect to find pathways related to diabetes and obesity diseases. However, I may also find in multiple patients that some pathways are recurring in the analysis. This would suggest a link to theorize about.

Yes, fair enough. I guess you can skip controls and assume simple over-representation. Just remember that in a WGS sample you expect ~5 million variants (including intronic and intergenic). Personally, I would use a collection of known healthy samples as controls, but you don't have to, you're right.

The concept of looking for enriched pathways makes sense, but you can't do it using the same tools you would use for gene expression data because the biases are different - longer genes are far more likely to contain variants by chance than short ones, and this needs to be accounted for in any enrichment analysis.

But a mutation can also be silent on the gene. It should also be taken into account, right?

Yes, although you'd probably want to take that into account as well.