text:
How to know if ID correspond to nucleotide sequence or protein with a script
I have a text file that contains a list of IDs (314 sequences):
AVP78031.1
AVP78042.1
ATO98108.1
ATO98120.1
ATO98132.1
...

My goal is to make a script (maybe using Python or Perl) to check in the list if all the IDs are nucleotide or protein sequences.
For example:
AVP78031.1 -> protein (this is a nucleotide sequence, I change nucleotide for protein to show an example).
AVP78042.1 -> nucleotide
ATO98108.1 -> nucleotide
ATO98120.1 -> nucleotide
ATO98132.1 -> nucleotide

Any idea to do a script?
Thank everybody!

If these are all GenBank or RefSeq accessions, you can use Entrez Direct for this as shown below:
$ cat accs.txt
ATO98108.1
ATO98120.1
ATO98132.1
AVP78031.1
AVP78042.1
$ cat accs.txt | epost -db nuccore | efetch -format acc
## no output because none of them are nucleotide accessions
$ cat accs.txt | epost -db protein -format acc | efetch -format acc
AVP78042.1
AVP78031.1
ATO98132.1
ATO98120.1
ATO98108.1

NOTE: This will work only if the accessions are currently live because epost does not find any suppressed accessions. For example:
$ cat accs.txt
NM_002826.3
NM_002826.4
NM_002826.5
$ cat accs.txt | epost -db nuccore -format acc | efetch -format acc
NM_002826.5

Here, all three accessions are valid nucleotide accessions but only the last one, NM_002826.5, is alive.
An alternate way is to use the accession prefixes defined here and come up with an appropriate regular expression query.

The alternate way (accession prefixes) is a great idea!


text:
Are there any databases for gene co-expression or expression pattern clustering?
I am currently working on gene clustering based on co-expression pattern in mouse brain. The problem is I do not have some solid way to test my result. Are there any suggestions for databases indicating groups of genes that co-express (in mouse brain), or even better indicating some similarity metrics between each pair of genes based on expression pattern?
Also, I do want some databases for differentially-expressed genes ;)
PS: By "differentially-expressed" I am referring to a pair of genes whose expression patterns in a single individual differ.

How about:

http://www.informatics.jax.org/expression.shtml
http://www.emouseatlas.org/emage/home.php

Other potential databases of interest:

https://biokeanos.com/search?q=co-expression+%2Bmouse

Thanks for your great answer!  The first two seems to be developmental expressions, while the third database is very helpful!  This site looks quite good: https://gemma.msl.ubc.ca/

source: Stack Exchange

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

sort: votes


text:
How to identify the DNA sequence of a gene in the complete genome sequence FASTA format file?
I need to identify the sequence of a gene in the complete genome sequence . I thought it was simple, instead it is not a straightforward task !
My method was the following: I downloaded the FASTA format file of the complete genome sequence of the archaea organism that I am studying, that is CP003685, and then I downloaded the CDS file from NCBI by clicking on ">Send to>Coding Sequencing>FASTA Nucleotide>Create File". In this way I have a txt file (always FASTA format file) in which I find a sequence for each gene (if I am not wrong).
So for example, I should find, in the complete genome sequence, the sequences of the CDS file (the last that I downloaded), right?
For example, if I consider:
>lcl|CP003685.1_cds_AFN03983.1_1007 [locus_tag=PFC_05180] [protein=DNA protection protein DPS] [protein_id=AFN03983.1] [location=complement(940023..940580)] [gbkey=CDS]
ATGCCAGAGCATAATAGGAGATTAGTTGAAAGGACAGGAATAGATGTAGAGAAGCTGTTGGAACTTCTCA
TTAAAGCGGCTGCAGCAGAGTTCACAACGTATTATTACTACACTATACTTAGAAACCATGCAACAGGTCT
GGAAGGAGAAGCAATTAAGGAGATTATCGAAGATGCTAGACTTGAGGACAGAAACCATTTTGAAGCCCTC
GTGCCAAGAATATACGAGCTTGGAGGAGAATTACCAAGGGATATTCGGGAATTTGCAGATCTAGCTTCTT
GTAGAGATGCCTACCTACCAGAAGAGCCCACTATAGAGAACATTCTAAAAGTCCTACTCGAAGCTGAGAG
ATGTGCAGTTGGAGTCTATACAGAAATATGCAACTACACTTTTGGAAAGGATCCAAGGACTTATGATTTG
GCCCTTGCCATATTACATGAAGAGATAGAGCATGAAGCATGGTTTGAAGAGCTCTTGACTGGAAAGCCGA
GTGGTCACTTTAGAAGAGGAAAGCCTGGGGAAAGTCCATATGTCTCAAAGTTTTTGAAAACTAGATAG

I should find in the complete genome sequence, this sequence, right ? But this is not. Why?
I suspect that I miss very important passages.
Thank you in advance.

See How to identify genes from a genome assembly of C. Elegans?. Exonerate and wise2 should be able to do this.

I simply copy the sequence of the gene (like the one I have written in the post) and in the file with the complete genome sequence (CP003685) I click ctrl+F and paste the reverse complement of the sequence of the gene and search it. But it does not find it. Sorry if I did not really specified this very crude way to do it :') . For the introns I have no idea...

It is a non-professional way but I do not understand why it does not work since the two sequences should be the same... the two files I am considering (as said below) are only two different ways to represent the genome DNA sequence, right ?

Well yes, kinda, but you have line breaks. You will never find the entire thing unless you convert it to a one-line format first. Try searching for a smaller section, say ~10 nucleotides and if the first one doesn't work try some others. You need to find a sequence that appears in a single line otherwise you won't find a match.

Thank you very much @terdon. I did not know this and now it seems to be finally clear putting all pieces together !

I would just blast it. When blasting locally, you need to first make a database from your genome, so assuming you got the command-line version of blast installed you can do something like
makeblastdb -in my_study_genome.fa -dbtype nucl
blastn -max_target_seqs 10 -db my_study_genome.fa -evalue 1e-10 -outfmt 6 -query my_downloaded_gene_of_interest.fasta -out gene_in_my_study_genome.blast

blast search for anything that is "similar enough" and also considers both strands, so for example, if there is an SNP in the middle of the gene. Hence even if the searched gene and the reference genome are from different sequencing runs / individuals, blast will still find it. Blast will also handle linebreaks for you (in fasta files the nucleotide sequence is quite often wrapped by 60 or 100 characters a line). Long story short, there is not much that can go wrong with blast (while with full text search you can easily shoot yourself in foot).
-- edit 2 --
As @terdon suggested, potentially if you would locating a transcript sequence (excluding introns), you would do better with software like exonerate or genewise, because they model exon/intron boundaries. However, even in those cases blast will give you an approximate position of the gene.

@Manuela this really should be done using a tool that models genes (try exonerate or genewise). Blast is great but will not be able to model introns correctly if your species has splicing.

In your sequence header you will see:
[location=complement(940023..940580)]

It's a reverse strand gene. Reverse complement the sequence, then search the genome sequence.

Thank you very much for the answer but it is still not working . I do not understand why. I used this software to make the reverse complement, copied and searched in the complete genome sequence but I still do not find it.

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


text:
How to understand methylation level?
For Bisulfite Sequencing data, we can calculate methylation level for a specific locus,
$$
\mbox{methylation level} = \frac{\mbox{# methelyated reads}}{\mbox{# total reads}}
$$
Here I think, for this locus, it can only have one status, methelyated or not. Or to say, methylation level here can only be 0 or 1. Why does it have other values in (0,1)?

You are right that a single molecule in a single position is either methylated or not methylated. However:
First, assuming your organism of interest is diploid (or of higher ploidy) one of the chromosomes could be methylated, the other not. That would give you a level of 0.5 and can be found in imprinted regions (where the paternally inherited chromosome is methylated and maternal not, or vice versa).
Second, except if you use a pure cell population of a single cell type it is likely that there will be cellular heterogeneity in your sample, e.g. from different types of blood cells or neurons vs glial cells in the brain. Some of these may have methylated CpGs in this locus (depending on cell-type-specific gene activity) while others can be unmethylated. Depending on the relative proportions of these cell types the methylation level can be anything between 0 and 1 for a mixture of cell types.
Thirdly, sometimes CpGs get grouped together. For example, 5 nearby CpGs are combined to give a single methylation level for this locus, and then the level can be 0; 0.2; 0.4; 0.6; 0.8 and 1.
Lastly, there could be technical artifacts or technical reasons why a position was not called as methylated (insufficient quality for a single read?), leading to deviation from 0 or 1.

Thank you very much! I noticed some databases supply a single methylation level for each promoter region, which I could not completely understand before. Now I think that's the third case in your answer. Your answer is really helpful.

Thanks for the explaination. One further question though, why would the methylation distribution are biomodal when 1) methylation could be heterozygous and 2) we have mixture of cell type where each type represent different methylation pattern?


text:
Bootstrapping cell type identifications (scRNA-seq)
I have scRNA-seq data held in a SEURAT object, and cell-types were identified with SingleR (essentially a vector of strings - Monocytes, Granulocytes etc.). I have two conditions - treated and untreated. I want to know if the difference in the relative frequency (%) of a given cell type between the conditions is significant. For this, I thought to perform bootstrapping, to obtain frequencies for many different 'simulated' datasets - and then perform a statistical test. However, I cannot understand how best to do this. Should I subsample the scRNA-seq data, with replacement and calculate the frequencies each time? Can somebody provide a step by step workflow on how to do this?

I would use permutation:
i) Get the difference that you are interested in, i.e. number/proportion of monocytes in one group vs the other.
ii) Get the relevant columns (condition, cell label, ...) from the @meta.data slot of your seurat objects and rbind().
iii) Shuffle the cell labels randomly and split the data frame into two (per condition).
iv) Compare the difference between the cell counts/proportions of interest from the two permuted data frames and save this value.
iv) Repeat iii and iv for 100s of times and calculate how many times have you observed a difference at least as extreme as your real observation (i). Once divided by the number of trials, this would give you a p-value stating how likely it is that you get the observed value.
I have to add that batch effects are severe in scRNA-seq and hence a linear model where you might be able to account for batches could be more helpful.


text:
Installation of PRANK MSA in WLS Ubuntu 20.04 LTS
I want to install PRANK on the Windows 10 Linux subsystem (Ubuntu 20.04 LTS), I have followed the installation instructions to no avail.
wget http://wasabiapp.org/download/prank/prank.source.140603.tgz
tar xvzf prank.source.140603.tgz
cd prank-msa/src
make
./prank

[ cp prank ~/bin ]
[ sudo cp prank /usr/bin ]

I followed those instructions but got stuck with cd prank-msa/src I don't know why they want me to enter a directory that doesn't exist. I also followed this video on how to install linux files but it wasn't helpful. Any suggestions would be appreciated.

Their instructions are outdated and don't match what exists on their server. The problem you're running into is not at the cd step but even before that - at your first step even. Your wget fails because that file does not exist. Go to http://wasabiapp.org/download/ and see for yourself - there is no /prank/ directory in there.
I tried using curl, and it wrote the 404 page as a .tar.gz file, so if you did that, the tar step would error out.


text:
Within and between sample count normalization
recently I came across a situation where RNASeq sample quality was expresses trough DESeq sizeFactor. So authors reported quality of their samples with respect to some publicly available datasets as good/bad if the scaling factor of their samples was equal/higher/lower than 1. What is more, they computed it using FPKM's from Cufflinks and executed geometric mean calculation over genes (as the one done in DESeq after log transforming the raw counts). Underlying logic being that with FPKM the samples were normalized "within" and with the additional DESeq, normalization "between", thus ready ready for any downstream analysis (differential expression, AS, etc.)
I haven't seen something like that before nor the use of DESeq scaling factor for estimating the quality of a sample (all samples are from the same species/tissue/condition/etc. )
Can anyone comment on why is this ok/not ok because I cannot wrap my head around the logic (then again I am just trying to replicate the study, no value contribution on my end.)
Thank you

could you please link to the paper in question?

The approach you are describing seems very strange.
Crucially, the Vignette for DESeq2 states that the model only works correctly with unnormalized counts as input:

It is important to provide count matrices as input for DESeq2’s
statistical model (Love, Huber, and Anders 2014) to hold, as only the
count values allow assessing the measurement precision correctly. The
DESeq2 model internally corrects for library size, so transformed or
normalized values such as counts scaled by library size should not be
used as input.

Normalized counts like FPKM are not supported and should be converted to counts using the appropriate import funcion.
I don't really understand the use of the normalization factors for quality control because the number of reads in a sample does not correlate well with the quality. Maybe you could link to the paper in case there is something hidden in the experimental design that would explain this.
If you are interested in quality control a useful metric in my experience is the 3' sequencing bias that can be calculated using, for example, RNA-SeQC.

The only reasonable use of scaling factors in QC is to ensure that they're not too out there, since if they differ by around 10x or more you likely had a sample that failed (you would pick this up in a normal PCA though) and need to exclude a sample.

Agreed, but a simple look at the number of reads per sample would tell you the same.

It sometimes doesn't, since it's often the case that you still have a number of highly expressed genes but little else. But yes, often this will be obvious when looking at MultiQC output or similar.


text:
PDB residue numbering
The residue numbering in PDB (in my understanding) coincides with UniProtKB align. However, this might not be consistent with that in PDB fasta file. To complicate it further, UniProtKB align can be a few mutations away from the PDB sequence. Disordered/unidentified residues in structures also create gaps in numbering.
Is there a way/suggestion to create a matched residue numbering between sequence and structure? For example
12345678 (residue numbering)
ATPWQMSQ (sequence)
---W-MSQ (PDB residues) "-" for missing residue in structure

Unfortunately, PDB numbering can be one of different schemes. The main three:

Uniprot numbering (canonical isoform)
Construct numbering with or without tag removal
Start from 1 for what has density

However, the are many many corner cases. In Proteopedia there is a lovely article about weird PDB numbering.
The best way to find the details about what part of the sequence is covered is to use the Swissmodel API. It's actually better than either the PDBe and RCSB PDB APIs. Swissmodel is a homology modelling service that maintains close threaded models (>30% ID), but also has PDB data. Here is a Python3 method I use, (provider='pdb') is what is required here).
def get_data(uniprot: str, provider='swissmodel') -> List[dict]:
        """
        Gets data from Swissmodel. Convert each entry with convert_SM_entry
        """
    assert provider in ('swissmodel', 'pdb'), 'Provider has to be swissmodel or pdb'
        url = f'https://swissmodel.expasy.org/repository/uniprot/{uniprot.strip()}.json?provider={provider.strip()}'
        data = requests.get(url).json()
        return data['result']['structures']  # [0] # from to coordinates description

The returned dictionary contains a start and end of the uniprot <--> PDB. It does not however, say where the missing density residues are. But in deposited structures the PDB numbering skips where these residues ought to be... So might not be an important problem.
Another database of interest is SIFTS.
To make a multiple sequence alignment of a non-deposited PDB, you need to extract the sequence and do the alignment here is a snippet using pymol module (from conda):
import pymol2

pdbfile = '' #your pdb file name..
selector = 'chain A' # or whatever chain you want
with pymol2.PyMOL() as pymol:
    pymol.cmd.load(pdbfile)
    fasta = pymol.cmd.get_fastastr(selector)

seq = ''.join(fasta.replace('?','').split('\n')[1:]) # remove header and non AA residues

Biopython has a dozen sequence alignment tools. But the original one still works fine
from Bio import pairwise2
from Bio.Seq import Seq

alignment = pairwise2.align.globalxx(''.join(seq).replace('?',''), "ACQQG")[0]

alignments is a tuple of the MSA sequences.
If you want to do the latter without coding. Extract the sequence from the PDB using PyMOL the normal application and load the pdb and then print cmd.get_fastastr('chain A') in the command line part of the GUI. Copy this. Go to Muscle or other onliner aligner and align!

Thanks @MatteoFerla. The swissmodel API is an interesting one. However, since the rcsb fasta does not have a UniProtKB id, perhaps I have to do an additional alignment to obtain the UniProtKB id? I am surprised rcsb has the sequence data (with the alignment, UnitProtKB id, missing residues, etc), but there is no API/downloadable available.

To go from PDB 4-letter code to Uniprot accession use the SIFTS API or mapping table from their FTP site. Or look manually at the PDB entry.
Do note that a uniprot record has an accession number, say for HRas P01112, and a entry name RASH_HUMAN. The "id" word is an ambiguous and can refer to either, but use to mean the latter. To confuse further, the entry name can change if it is moved from Trembl to Swissprot —say the genome of that organism was sequenced and before the protein was known from an EST read or something.

The RCSB PDB has changed recently, dropping a few annotation from its protein viewer, which I am guessing will be re-added —once that happens I assume a sensible API will appear. How it used to work was that HTML/JS chunks were assembled in the backend, so there was no API for uniprot data and you could not reverse-engineer an ajax call into an API.

Thank you so much. It is indeed a very comprehensive explanation.


text:
RCSB API: limited returned result
Is there a limit of query results on RCSB API?
Interestingly, no matter how to adjust my query criteria, only 10 assemblies are returned. Here is an example query json: (truncated from a longer query file)
{
  "query": {
    "type": "group",
    "logical_operator": "and",
    "nodes": [
      {
        "type": "terminal",
        "service": "text",
        "parameters": {
          "operator": "greater",
          "value": 1,
          "attribute": "rcsb_assembly_info.polymer_entity_instance_count"
        }
      }
    ]
  },
  "return_type": "assembly"
}

Example python code requesting query:
url_query = 'https://search.rcsb.org/rcsbsearch/v1/query?json='
query_str = open(query_json, 'r').read()
query_str = query_str.replace('\n', '').replace(' ', '')
url = url_query + query_str

r = requests.get(url)
content = r.content.decode('utf-8')
result = json.loads(content)
print(len(result['result_set']) # always equal to 10

Maybe ask RCSB at https://www.rcsb.org/pages/contactus.

Staff from rcsb kindly reminds me the solution - add the following lines in the json query.
"request_options": {
"return_all_hits": true
}


text:
DNA seq cutting using Linux
I need to use the cut command on a human insulin sequence row wise but i don't know the command.I already ran byte wise command but the gave me column wise results. I cant use the -f option since i need to use a delimiter with it and the seq does not have any delimiters.
please answer asap

Could you please show the commands you tried so far and how the result should looks like? If you could give a simple example of how your data, it would be a good help to understand your problem.

cut is a linux command, but might not do what you intend. Could you please describe your problem with greater clarity?

I have an assignment where i have to use Linux commands on a bioinformatics problem i choose  to cut sort and uniq the sequence. however the cut command worked for byte wise cutting but not field wise since fields dont exist in a seq. also can someone tell me why do we sort a seq because wouldn't that result in a mutation.

Use either generic *NIX cut (you do not need a delimiter if you use -c option) or head tools. Even better, use any of a number of bioinformatics tools for sequence manipulation, such as EMBOSS extractseq.
EXAMPLES:
# Write into file output.txt characters 1-20 from each line of file input.txt:
cut -c1-20 input.txt > output.txt

# Write the first 5 lines of input.txt:
head -n5 input.txt > output.txt

# It is easier if you have your sequence in fasta format.
# Write characters 1..20 and 100-110 from input.fasta to output.fasta:
extractseq input.fasta output.fasta -regions "1-20 100-110"

One of the easy ways to install EMBOSS suite is using conda (from bioconda channel).
SEE ALSO:
EMBOSS extractseq help: https://www.bioinformatics.nl/cgi-bin/emboss/help/extractseq


text:
Blastp MSA to the same length
To generate aligned protein sequence truncate to the same length, I use blastp in biopython for MSA. I have multiple queries and one subject sequence. However, the alignment never returns the same length.
How can I align all sequences altogether, thus resulting in a set of sequences of the same length?

Please consider accepting/upvoting the answers below.

@Matteo Ferla wrote in the comments,

It sounds like you'd want a MSA of the whole set and not a set of pairwise ones. If global MSA is what you actually want, try ebi.ac.uk/Tools/msa/muscle. Then to trim the ends before phylogeny there are dozens of tools.

The reason is because Blast aligns locally, i.e. a window, while a global alignment is the whole entirety of the sequence

BLAST isn't appropriate for multiple sequence alignment.
BLAST generates pairwise alignments. Furthermore, the aligning regions shown in a BLAST report only include subsets of the hit sequences.
To create multiple sequence alignment of hits you want to:

run blast
extract the hit sequences
perform multiple-sequence-alignment (e.g., using MAFFT)

Extracting hit sequences can be done using the command-line. Check the docs for blastdbcmd.
If you prefer a point-and-click interface, you can also do this using my lab's SequenceServer BLAST software by simply selecting the hits that you want to retain, and in the Download section, clicking "FASTA of selected hits". Or if you want all of them, click "FASTA of all hits".

Please add a statement stating that your lab created this software that you're posting about on multiple threads.


text:
Parse Error while using AutoDock Vina
I am trying to dock ligands with proteins to estimate the docking potential however, I am facing this error while docking. Any help will be highly appreciated.

could you share the content of the line 574 of the Spike_C.pdbqt file? You can use sed for select one line after and one before of the 574 specific line  sed -n '573,575p' Spike_C.pdbqt

Could you also let us know how you converted your ligand to the PDBQT format (prepare_ligand4.py, obabel etc.)

Also (minor point), having the grid box the size of your protein (270_000 Å) is not a good idea: you may know where the ligand binds —i.e. receptor binding domain. In the case of a VLS you want scores for it to bind there. Otherwise a lot will stick to the energetically unhappy transmembrane domain doing nothing.


text:
Error while splitting PDBQT file
I am trying to split to a pdbqt file containing 178 compounds into individual pdbqt files. I have been using the Autodock Vina_split.exe file to split the multiple entry file however I keep running into the following error. What can I do to prevent the error.

(A) Empty file? Does the file open in PyMol? Also, you can always split it there with split_states.
(B) Is the file really called output.pdbqt or is it output.pdbqt.txt and Windows is hiding the extension?

There could be many reasons. Maybe lack of permissions or the space in the path. Can you open the file with a different program? Are you on Linux or Windows etc.


text:
How do I create Protein Ligand images
I was wondering how I can create images showing the bonds and the residue codes involved in the bound Protein Ligand complex obtained after Molecular Docking as depicted in the image below.

Pymol, in either its commercial or open source edition, will be sufficient to make figures such as your examples. Pymol recognises "polar contacts" for the intermolecular interactions and has labels that can be used for residue codes.
It is also scriptable and is widely used.

thank you very much for getting back to me and helping me out.

Yes, PyMol is one of the best programs for rendering IMO. Usually, you would create folder with your PDB file and a file that contains all the commands to create your image like shown here. You can create more than one file. Because PyMol is based on Python you can use all python syntax here. E.g. you can create multiple images and with for loops you can create series of images that you can combine to a movie later on.
Other programs to check out would be VMD. I use VMD for quick videos and images. The PyMol images I like better most of the time.
Or newer software like Samson modeling platform. I am not sure about the licensing conditions. Have not tried it yet.
Biopython does not support the generation of 3D molecular images. It is more a library for sequence data.

Great answer. On the biopython point, NGLView is Jupyter notebook visualiser which can accept an impressive multitude of structures, including Biopython structures —however, this is totally not suggested for paper figures (as the OP wants).

@soren thank you very much for getting back to me and helping me out.

The other contributors have given great suggestions. I thought I ought to mention Michelaɴɢʟo, which is an online PDB file/PDB code/PyMOL-PSE/gene-name visualiser, where one can create an interactive page, add a description with special protein links (prolinks) that change the interactive protein view and then share the URL securely.
It has a modal to make the selection easier, but uses the NGL selection language and to show a ligand, one would add a residue type prolink selecting ligand. If the neighbourhood-showing radius setting is greater than 1.0 Å (default is 3.0 Å) one will likely see hydrogen bonds with the neighbours.
One can try it without registering (free), but cannot keep the page for longer than a month or edit the description.
Disclaimer: I am the author.

Thank you very much for getting back to me and helping me out.


text:
Hybrid approach in metabolomics (NMR and MS)
Is it currently common in metabolomics to combine two quantification methods, namely NMR (nuclear magnetic resonance) and MS (mass spectrometry, either liquid or gas)? If so, what would be profits from such an approach? Is it a promising direction or people are going to stay with a simple approach and they will be sticking to either of these two?

Hi and welcome, perhaps you could explain whats yuor investigation?

Hi Michael, I am just curious if running an NMR experiment will boost the insights that I may obtain from pure MS experiment. I have the feasibility to run both of them, I am just wondering if it makes sense.

It strikes me as more of a wet-lab question than bioinformatics, but its just my opinion

Not necessarily, wet-lab question would be, for instance, about applicability of a given method to the samples that I have. But I know the pros and cons of both methods. What is of my concern is the next step - whether the analysis utilizing both of these brings something more to the table in comparison with sole MS / NMR.


text:
Get gene sequence based on the annotation
I've got the reference genome with Python like so:
sequences_by_chr = {}

with gzip.open("data/Zea_mays.B73_RefGen_v4.dna.toplevel.fa.gz", "rt") as f:
    for seq in SeqIO.parse(f, "fasta"):
        sequences_by_chr[seq.id] = seq.seq

I've also parsed a GFF3 annotation and would like to get the gene sequence based on said annotation.
Suppose the following:
chromosome = 1
strand = "-"
start = 45286331
end = 45290076

Since it's on the "minus" strand, should I do the following:
sequences_by_chr[chromosome].complement()[start:end]
or should I use the reverse_complement()?
I'm completely confused by this, would appreciate any hint!

You can parse the GFF3 to generate a bed-like file with gene ranges and use that with bedtools getfasta to obtain a multi-FASTA file with gene sequences. An example case is shown below:
## get fasta and gff3 files
wget ftp://ftp.ensembl.org//pub/current_gff3/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz
wget ftp://ftp.ensembl.org:/pub/current_fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## unzip fasta
gunzip Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## generate genes.bed
zgrep -v '^#' Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz | awk 'BEGIN{FS="\t";OFS="\t"}($3=="gene"){print $1,$4-1,$5,"name",1000,$7}' > genes.bed

## create gene fasta
bedtools getfasta -fi Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed genes.bed -s -fo genes.fa
```

Thanks, I'll keep that in mind! But at the moment, I'm operating on anair-gapped server where I don't have the tools and can only use Python.

In that case, you may want to check out pybedtools or pyfaidx. In both cases, you should be able to provide a list of ranges along with an indexed FASTA to extract sequences in multi-fasta format.

You should do:
if sequences_by_chrom[strand] == "-":
    rev_comp( sequences_by_chr[chromosome][start:end] )
else:
    sequences_by_chr[chromosome][start:end]

where
def rev_comp(sequence):
    d = {'A':'T', 'T':'A','C':'G','G':'C'}
    return ''.join([d[i] for i in sequence][::-1])

DNA sequence is always provided on the plus strand, if the gff tells you it is on the minus strand, you need to reverse complement it, because the start of the gene in your example will be 45290076 and the end will now be 45286331.
I can give you an example below, this is a potential coding sequence from yeast, no introns and on the minus strand:
gunzip -c test.fa.gz
>1
TCATCCCCTTATATAAGATGGGAAACGAAGGTAAAAGATTAACGATAGCAT

If we read it in like you did, and if you apply the complement, this is not a coding sequence.
chromosome="1"
start = 1
end =51
sequences_by_chr[chromosome].complement()[(start-1):end]

Seq('AGTAGGGGAATATATTCTACCCTTTGCTTCCATTTTCTAATTGCTATCGTA', SingleLetterAlphabet())

Getting the sequence followed by reverse complement gives you a coding sequence, starting with ATG ending with TGA:
sequences_by_chr["1"][(start-1):end].reverse_complement()

Seq('ATGCTATCGTTAATCTTTTACCTTCGTTTCCCATCTTATATAAGGGGATGA', SingleLetterAlphabet())


text:
Need a alternative or more complex version of venn diagram in python for matching dna sequences
I am in google colab and I have combined and set up a data frame list of sequences that goes like this
Location ID Sequence
1.1     ........ A ........ AAGAGATA
1.2     ........ A...........   CGGATGAT
1.3     ........ B................   ACGTGATC
1.1     ........ B ........ AAGAGATA
1.2     ........ B...........   CGGATGAT
1.3     ........ C................   ACGTGATC
So I need to make a alternative to Venn diagram because there are about 20+ ID types within a couple hundred rows of list just like this. And some of the sequences repeat or are unique to its id so I need something to be like
CGGATGAT Shows up 4 times in B,C,A, and D
ACGTGATC shows up 1 time in B
AGTGTCTC shows up 5 times in E,A,B,C, and D.
And then have some kind of graphical visualation to show all of these common sequences and where they are located and which are unique. Using python. This might be asking for too much but I'm not sure. I would use something like vennr to make a bigger complex version of this but I don't really know R.

Something like this, perhaps? https://64.media.tumblr.com/tumblr_m723nmp5OM1qbh26io1_1280.jpg (http://www.promusa.org/blogpost174-The-best-genomics-Venn-diagram-ever-deconstructed). Bear in mind that there is likely no way to visualize this that will be as effective as a table of numbers.

To get around limitations in using Venn diagrams for set overlaps, when there are more than three sets, back around 2013 or so I created something I called an Eulergrid plot (example at the bottom of the page), which an UpsetR plot appears to recapitulate, today.
The Eulergrid code I wrote was in a mix of Perl and R; the UpsetR plot code uses R. There appears to be a Python port of UpsetR that might work, as well.

The Python upsetplot module works pretty well!


text:
Produce .vcf file of ALL mutations in .bam file
I got tasked with a problem, that i thought would be quite simple to solve, but turned out to be quite tricky. Our lab is running targeted mutagenesis experiments in yeast using crispr base editors. What we have done is set up an experiment where the base editor selectively mutates a region inside of a gene of interest in a culture of yeast cells. We then extract the DNA of the whole culture and perform sequencing of our gene of interest (the size of which is around 1kb).
What i would like to do is plot the mutation frequency of say, C -> G edits, across the length of the entire gene (around 300bp), with the hope that i see a spike in mutations at the site where the base editor binds. Or at least, higher mutation rates at this site compared to the background mutation rate.
These edits are quite rare (we think), and since it's in a culture of multiple cells, they don't occur at the exact same position on the gene.
I initially thought i could do this by quality filtering the reads, aligning them to the gene sequence in order to produce a .bam file, and then piping this file through variant calling algorithms such as vcftools. However, i noticed that i get very few variants out of this, even when using extremely relaxed settings (e.g. p-value < 1 for the variant call).
Im therefore wondering if anyone knows how to extract from a .bam file ALL the mutations across all reads that do not match the reference sequence?
I'm aware that this will probably lead to quite a huge .vcf file, but since the reference sequence is only 300bp, it should still be manageable? It could also be the the experiment does not work, but i would at least expect some mutations in the sequencing reads compared to the reference gene sequence.
Any help would be greatly appreciated!

Ensure that the variant caller stringency is relaxed all the way. For example in freebayes use --min-alternate-count 1 --min-alternate-fraction 0 .

You could perform an mpileup with samtools. I think this would be the most sensitive approach. Here's the man page: http://www.htslib.org/doc/samtools-mpileup.html
You might want to disable BAQ correction. If you specify a reference sequence, then any letter in the fifth column would represent a non reference base. Reference bases are coded as . or ,.
The reason why vcftools doesn't work is probably because it assumes a diploid individual, which is not what you have. You have lots of samples mixed together. Freebayes allows relaxation of the diploid assumption. Here is the relevant example from the documentation https://github.com/ekg/freebayes :

Generate frequency-based calls for all variants passing input thresholds. You'd do this in the case that you didn't know the number of samples in the pool.

freebayes -f ref.fa -F 0.01 -C 1 --pooled-continuous aln.bam >var.vcf

You could try a cancer/normal variant caller like mutect, but I don't think you'll have much luck.
The most straightforward approach is to use samtools mpileup and to process the output yourself. You could also use bcftools mpileup to get a VCF, and then extract the allele count from the INFO field. If you are a Python lover, then you could also use pysam to both perform the pileup and process the results.


text:
How to find possible viral subsequences in a Eukaryote genome?
I have an assembly of an algae and want to find if it has anything analogous to endogenous retroviruses' subsequences in its DNA.
My guess is I should manually blastn all subsequences until I find (or not something) that could have come from a virus.
Is there a better way to do this?

Is it your analysis any different to the detection of horizontal gene transfer? (Or at least I would search for this keywords for an answer)

you can try with kraken: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-3-r46

Your guess sounds reasonable to me but you should consider using KRAKEN2. You can create your custom database (you would have done it with makeblastdb anyway) and kraken will check k-mers against the database to find homologous sequences.

I don't feel Kraken is good for this task. Kraken taxonomically classifies each sequence, while the signature of integrated viral DNA will be a single sequence with both host and virus parts

@Chris_Rands. Preparing your custom dataset with only viral genome and with the right k-length for the kmer analysis is the only method I know of existance today. You can use manually blastn (honestly, I do this myself if it's not for publication) or let the job to a well developed (and fast) software out-there. Please, let me know if you think of a better or potentially better solution!

Your suggestion might be fine for detecting viruses in general, but not integrated viruses specifically- this requires a fundamentally different approach because you are looking for chimeric reads

@Chris_Rands chimeric reads could be found using shorter kmers (at the expense of dropping specificity). Do you know of any better or alternative solution?

Searching directly for a reverse transcriptase may find recently incorporated retroviruses.


text:
How to identify to each scaffold a read belongs to, inside a .sam file?
I have a fasta file assembly and combining it with the raw reads we produced a .bam file which I converted to .sam .
The .sam information lines look like this:
A00321:42:HLLVYDSXX:2:2302:6153:3505    99      NODE_1_length_3415511_cov_137.721502    16      60      128M    =       607     742     CGATTAGTCCGGCCAAATCGCCGTCGAGCGCAATGAACATAACGGTCTTGCCCTCAGCGCGCAGCGCATCGGCCTTGGCGTCGATTGTGGAGTGCTCGACGCCCATGATGTCCATCATAGCACCATTG        FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        RX:Z:TTGAGGGTATAGTAGT   QX:Z:FFFFFFFFFFFFFFFF   TR:Z:GACACCG    TQ:Z:FFFFFFF    BC:Z:AGTTGCAG   QT:Z:FFFFFFFF   XS:i:-10        AS:i:0  XM:Z:0  AM:Z:0  XT:i:1  RG:Z:over_1kb:LibraryNotSpecified:1:unknown_fc:0        OM:i:60

Separated by mandatory fields it would be something like this:
QNAME: A00321:42:HLLVYDSXX:1:1644:2248:3881
FLAG: 99
RNAME: NODE_1_length_3415511_cov_137.721502
POS: 1
MAPQ: 60
CIGAR: 1S127M
RNEXT: =
PNEXT: 536
TLEN: 386
SEQ:  ATCGGGTCTGACACCGCGATTAGTCCGGCCAAATCGCCGTCGAGCGCAATGAACATAACGGTCTTGCCCTCAGCGCGCAGCGCATCGGCCTTGGCGTCGATTGTGGAGTGCTCGACGCCCATGATGTC
QUAL: FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

I'm actually interested in the meta data. I want to know how the RX: and BC: fields are distributed across the scaffolds in the original assembly.
I imagined the .sam file already contains the information about the assembly used to produce it. If I'm wrong, I'm sorry and please correct me, I'm just assuming.
What I want to do is, for each read in the .sam file, I find out its position in the assembled scaffold, and I record, Read_ID,Scaffold_ID,Read_Position_Inside_Scaffold,RX,BC
Then I want to use that database to analyse the distribution of RX and BC inside each scaffold.
That's what I want.
Ultimately what I'm trying to do is evaluate the quality of my assemblies based on the Barcode distribution.
I'm good at programming and parsing, I'm just having trouble figuring out, where, inside the .sam file, can I find the scaffold and scaffold position of each read.

What I want to do is, for each read in the .sam file, I find out its position in the assembled scaffold, and I record, Read_ID,Scaffold_ID,Read_Position_Inside_Scaffold,RX,BC

I'm just having trouble figuring out, where, inside the .sam file, can I find the scaffold and scaffold position of each read.

You already listed all those:

ReadID -> QNAME
ScaffoldID -> RNAME
Position -> Pos

Assembly evaluation is a wide topic, but other software you can use includes BUSCO, QUAST or something like LTRi. Here's a more in-depth guide by the SciLife lab in Sweden. Usually, when aligning reads back to an assembly, you care about the mapped insert size distribution compared to the theoretical insert size distribution, which tells you about collapsing or exapnding regions (i.e. structural misassemblies). This is something that tools like FRC_Align do.


text:
In Galaxy, using the bracken tool, how do I get a Kmer distribution?
As you can see in the image, the Kmer distribution field is empty: "No options available"
What am I missing?
Is there another tool I must use first to generate a Kmer distribution?

Try posting here instead: https://help.galaxyproject.org/


text:
What to use to determine SNPs in a population? How to find significant SNP variation between generations?
I have a population in the hundreds and did an Illumina pool-seq (ended up with ~100bp reads for the whole population).
I have the reference genome for the species and now need to align all the population reads to the reference genome, figure out which are SNPs and finally determine which ones are changing from one generation to the next (about 20generations difference).
Right now I'm trying to use popoolation2 https://sourceforge.net/p/popoolation2/wiki/Main/
but it seems to be a little old as it has not been updated in ~10 years.
In Galaxy, I couldn't find any pipelines dedicated to determine population SNP variation through generation.
So two questions:

what software to use to determine significant variation in SNPs in a population (in poolSeq)?
can you suggest me a pipeline for this? if one exists.

First focus on calling SNPs accurately. See the BWA, SAMtools, BCFtools pipeline (as described here) or GATK4 (as described here).
Once you have the [BV]CFs, then just code up the analysis to see changes over generations.


text:
Simulation software recommendation to confirm pipeline validity
I'm trying to detect convergent evolution in different samples of the same species.
I'm just trying to see if past 20 generations the nucleotide frequency of snps under pressure converges somewhat for all samples.
What i want now is to validate my pipeline using a simulation.
I would like to simulate the evolution of populations and see if the pipeline i built would detect the snps from the reads.
I'm using pool-seq and illumina with reads of around 100base pairs.
Can you recommend me some software and maybe point me to the tutorial?


text:
Building a group representation of genes' ontologies
I have a list of genes (about 2000) from lines of drosophila subobscura, in the following format:

LOC117900589
LOC11788959
LOC11790331

I'm trying to represent this set of genes in a useful way.
So, I was recommended to represent the ontology function of these genes in a graphical way.
I was recommended this tool:
http://cbl-gorilla.cs.technion.ac.il/
One inputs a list of genes and it outputs this kind of graph:

The problem is that, this tool, GOrilla only accepts genes from drosophila melanogaster.
I have to either:

a) convert the genes id to a homologous gene in melanogaster
b) find an alternative tool

A solution to either a) or b) would be helpful. Or an alternative solution.
For a) I tried:

Ncbi's homologene: which doesn't seem to be working for a the genes ID I tried so far (I didn't try them all)
Blasting the gene sequence on melanogaster's genome and looking through the list manually to find an homologous gene.

So far these methods have been impractical for the number of genes I have.

Have you tried DAVID for functional annotation (Gene Ontology and KEGG terms)?
I usually use DAVID to extract GO and KEGG terms and plot the results using R.
If you also want to plot a protein-protein interaction (PPI) network for your gene list you can use STRING database to extract interactions and network and use Cytoscape for visualizing this result.

I tried putting a bunch of entrez id genes and I got nothing useful it seems? https://imgur.com/tZc6Vk8.png


text:
How is an X chromosome encoded into a fasta string?
A human "X" chromosome has a centromere and two "identical" chromatids. If the chromatids are not identical, this fact is not assembled, correct?
The fasta string for a chromomsome is actually representing just one of the chromatids, correct?
So from my understanding, for diploids, the fasta chromosome is actually a mix of 4 chromatids: 2chromatids * 2 pairs of chromosomes into one string. Which for some heterozygous nucleotides (same percentage of each allele) will just be randomly selected from the reads as the correct one.
(unless we are talking about some methods like 10x that can distinguish between chromosomes pairs)
I was thinking which side of the X is read first into the fasta string, and I came to the above conclusion.

So, the two chromatids you see only during meta-phase of replication. In more "regular" state of a cell, there are only the two haplotypes (in the case of a diploid organism)

Well, and regardless of the ploidy, we usually aim for a haploid (monoploid) representation, so it is usually exactly as you say - an arbitrary one of the two states on each heterozygous site is chosen.
However, there are also approaches to phase genomes (to sort which alleles are on which haplotypes). That is possible not just with 10x, but also with long reads (to certain extend) or using trio-binning - a strategy to sequence two parents and the offspring.
However, even phased genomes still contain the reference (fasta) in a haploid state and usually provide the alternative haplotype as a second fasta file.
For the record, these apply to all chromosomes, not just the X.


text:
Tracing and quantifying inheritance over many generations
Apologies if this is the wrong location for this question, and please feel free to point me in the right direction with appropriate Anglo-Saxon embellishments...
I'm after some ball-park numbers for the following questions (plus any hints how to calculate them for different params - e.g. less generations, larger starting pop.):
Given a starting population of 10,000, and given 2,500 generations, with an end population of 5,000,000,000

What are the chances of being a direct descendant of 1 individual of the original 10,000?
If you had a DNA sample of that individual, would you be able to tell if you were a direct descendant?
Would siblings of the original individual confuse things with regard to Q2?

I don't know how to solve this, but your question reminded me of subjects discussed in the book 'The Ancestors' Tale' by Richard Dawkins. Chapter 0, which deals with human(oid) population genetics, is the relevant part. At least a good part of it appears to be available for free on GoogleBooks https://books.google.co.uk/books/about/The_Ancestor_s_Tale.html?id=s86_CQAAQBAJ. Maybe it'll help?

Thanks Laura. If I've read that right, then the answer to 1 is "assuming they have any living descendants, then it's 100%" (which seems reasonable).

Err, how do I flag a question as 'answered'?

As for question 2, you could probably infer whether or not they are a direct descendent by looking at the uniparental markers (either Y chromosome for men, or mtDNA for females)

Hi @Laura if you place your comment as an official "answer" Tom can then "accept it" and we can all upvote you.

Since it seems my comment helped, I've been suggested to post this as an answer...
I'd reccomend everyone interested in these sort of population genetics to have a look at 'The Ancestors Tale' by Richard Dawkins. In particular, chapter 0, which deals with human/hominid genetics. That should help answer questions. The relevant part of the book appears to be available online for free in GoogleBooks (https://books.google.co.uk/books/about/The_Ancestor_s_Tale.html?id=s86_CQAAQBAJ).
Below I'll note down the answers from the comments section, not because I wrote them but I figured it might help anyone who stumbled over this question afterwards:

The chances of being a direct descendant of 1 individual of the original 10,000 is 100%, assuming they have any descendants. - Tom Melly

If one had a DNA sample of that individual, one could use uniparental markers (either Y chromosome for men, or mtDNA for women) to infer whether they are a direct descendant. - 4galaxy7

(My guess is that it doesn't have to be uniparental markers. Other genes might work just as well. I'm thinking of Neanderthal/Denisovan work. But I might be wrong.)

Question 1
If I understand the question, then you are asking 'Given a starting population of 10 000 individuals and 2 500 generations, what is the probability that one of the 5 000 000 000 individuals descendend individuals is directly descended from one of the founders?'.
If we are assuming no migration in or out of the system, then the answer must be 1, since there are no ways to be descended from someone outside the system. Please let me know if I have misunderstood the question.
Question 2
For this, it's only possible to use uni-parental markers. This is because inferring exact genealogical relationships using autosomal markers beyond something like 3rd cousins is extremely difficult, owing to the stochastic way in which recombination breaks up segments. Beyond a few generations, only the uniparentally inherited mtDNA and non-recombining Y-chromosome can be informative about relatedness.
Uniparental markers don't recombine, so mutations aside, you expect to inherit an exact copy of your mothers mtDNA marker. If you know the number of generations which are likely to separate you and the putative ancestor, you can estimate the some kind of likelihood of being directly descended from the ancestor given the number of pairwise differences between the individual and the ancestor.
This is what they did when they wanted to identify the remains of Richard III. It's a good paper and easy to read / open access, so I'd recommend checking it out.
Question 3
Simply, no. The inference of being a direct descendent from one of the ancestors would be independent from any of the other potential descendants in the sample.

Many, many thanks - just sorry I can't flag more than one 'answer'.


text:
Why is my nextflow script executed sequentially?
I have the following Nextflow script:
nextflow.enable.dsl=2
include { check_folders; downloadClinvar; downloadReferenceGenome; downloadBedFiles; downloadGatkVcfFiles } from '/home/user/Desktop/pipeline/download_data.nf'
workflow {
    check_folders()
    downloadClinvar(check_folders.out)
    downloadReferenceGenome(check_folders.out)
    downloadBedFiles(check_folders.out)
    downloadGatkVcfFiles(check_folders.out)
}

The graph that is generated for me using the -preview -with-dag flags is the following:

When I run the workflow it works fine but it performs each step sequentially, and I want it to download everything at the same time.
What could be happening to me? Why is my nextflow script executed sequentially?

Cross-postd https://stackoverflow.com/questions/78445082/why-is-my-nextflow-script-executed-sequentially

Does the machine you run on have sufficient resources to run these processes in parallell?

Yes Pallie, I have Intel i8 7th generation, with 16GB of RAM

If your processes request a large amount of RAM/processors per process, that still could be the reason you're not getting parallellised processes


text:
Statistical or performance test for genomic intervals between 2 samples
Is there an appropriate statistical or performance test to perform when I want to compare genomic intervals between 2 samples. Sample A is the the benchmark of comparing the intervals and sample B is the sample in question. What I am really looking for is not only how many intervals from sample B did indeed overlap with sample A but additionally what intervals from sample B has intervals that don't overlap with sample A and also those intervals in sample A that don't overlap with B and being able to measure the performance of how accurate sample B intervals overlap with sample A. This in some sense, can be translated into the following:

True positive: intervals from sample B that overlap with sample A.
False positive: intervals from sample B that don't overlap with sample A.
False negative: intervals from sample A that don't overlap with sample B.

If I consider the above schema, then I will not be able to do any of the performance and statistical test (ROC curve or accuracy) because it is not possible to know what the true negative and the true negative are needed parameters. Then, what would should I be using to be able to incorporate the above schema and still apply performance and statistical testing.
I highly appreciate any help or pointers in advance.

Could you please give more context around your question? What samples? What gene? Can you provide any examples of what you're interested in?
Abstract questions about hypothetical situations are not easy to answer.

Are you looking for some coding solution to get which A interval fall into B interval ?

Please consider @MaximillianPress's for upvoting/accepting. Accepting an answer does help with site stats

I would recommend looking into BEDtools, assuming that you are talking about genomic intervals, which I think you are based on the question. Possibly a diagram would help.
BEDtools has a suite of statistical tools for evaluating overlaps, including Fisher Exact Tests and also approaches based on the relative distance of 2 sets of intervals. So as long as your genomic intervals are defined by BED files and a genome file, this should be quite easy, e.g.
bedtools fisher -a A.bed -b B.bed -g genome_file.genome
# will print p-values for left, right, and two-tailed tests


text:
Python equivalent to R GRanges and IRanges
Problem:
I am trying to convert some codes written in R to Python and part of that conversion process is find classes equivalent to the GRanges and IRanges from the GenomicRanges R package in Python. https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html#granges-genomic-ranges
I couldn't find any equivalent library written in Python that does the same operations as the aforementioned R classes. There is however pyranges library in Python but it is not as flexible as GRanges and doesn't have IRanges implemented in it. https://github.com/biocore-ntnu/pyranges
For example in R I could add the following in the GRanges class:
gr = GRanges(
seq = Rle(DF$chr),
ranges = IRanges(DF$Start,DF$End),
ref = DF$ref,
alt = DF$alt,
sampleID=DF$sampleID,
seqlengths = chromSize)

However, pyranges don't allow you to add arguments other than (chr, start, end, strand). However, in my case, I need to have the option to add sample_id, seqlengths and other arguments.
Thank you in advance.
Best

@TheUnfunCat, do comment or provide your take, since it's your package :)

I am not sure if your question has anything to do with IRanges. If I get it correct, the limitation you have with pyranges now is adding extra columns to the meta data, like adding values to @values slot in granges.
You can do it using insert or setattr :
import pyranges
import pandas as pd
import numpy as np
gr = pyranges.random(10)

setattr(gr,'alt',np.random.choice(['A','T','G','C'],10))

gr.insert(pd.DataFrame({'ref':np.random.choice(['A','T','G','C'],10)}))

+--------------+-----------+-----------+--------------+------------+------------+
| Chromosome   | Start     | End       | Strand       | alt        | ref        |
| (category)   | (int32)   | (int32)   | (category)   | (object)   | (object)   |
|--------------+-----------+-----------+--------------+------------+------------|
| chr1         | 27431195  | 27431295  | -            | C          | G          |
| chr2         | 197045893 | 197045993 | -            | A          | G          |
| chr4         | 86316012  | 86316112  | -            | C          | A          |
| chr8         | 4560598   | 4560698   | -            | G          | G          |
| ...          | ...       | ...       | ...          | ...        | ...        |
| chr11        | 130395526 | 130395626 | -            | A          | A          |
| chr12        | 28557156  | 28557256  | +            | T          | G          |
| chr13        | 55519337  | 55519437  | -            | T          | C          |
| chr14        | 1807028   | 1807128   | +            | G          | A          |
+--------------+-----------+-----------+--------------+------------+------------+
Stranded PyRanges object has 10 rows and 6 columns from 10 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.

Or directly adding them like this (thanks to @ChrisRands for pointing it out):
gr.alt = np.random.choice(['A','T','G','C'],10)

Thanks @Chris_Rands, I have added it to the answer

thanks alot for your answer @StupidWolf. What about mimicking the IRanges in R? does pyranges support that too?

what function do u need from IRanges? as in.. mimicking IRanges is really broad. Is there something you need to do specifically with IRanges and you want to reproduce that in python?

Mimicking this type of call in R, IRanges(start=1:10, end=11)?

whats the difference between that and creating a pyranges object with pyranges.from_dict({'Chromosome':1,'Start':np.arange(10),'End':11})

you cannot mimick a call... it's really weird... you have a function and it must do something, which is my question in the first place

Yeah, functionality can be "mimic"-ed i.e. with an equivalent implementation, but one cannot expect to Find-Replace function names and expect code to be ported over to a new language. I guess OP has some R code and wishes to implement that in py.

In pyranges you would not need a call like the one you show. You can just convert a dataframe with the columns you mention above to a pyranges like this:
import pyranges as pr
gr = pr.PyRanges(your_df)

This avoids the code duplication shown in your example. The column names in the df are preserved in the PyRanges.


text:
Error When Using biocLite as an installer in rpy2 python library
I am running into an error when I tried to run the following lines of code, the error at line 3:
1 from rpy2.robjects.packages import importr
2 base = importr('base')
3 base.source("http://www.bioconductor.org/biocLite.R")
4 biocinstaller = importr("BiocInstaller")
5 biocinstaller.biocLite("seqLogo")

rpy2.rinterface_lib.embedded.RRuntimeError: Error: With R version 3.5 or greater, install Bioconductor packages using BiocManager; see https://bioconductor.org/install
How can I use BiocManager installer instead of biocLite  and reflect that on my code? I couldn't find any documentation or source that talk about this error from within Python.
Thank you in advance

you need upgrade your R base version.

Definitely recommend installation via Conda as @zorbax points out. This normal stuff for installing rpy2 As a word of warning further errors with rpy2 are also normal and as a general rule of thumb rpy2 is a nightmare. For R's ggplot2 via rpy2 is outdated and will not import the full funcationality of the ggplot2 into Python3. This is irregular given have important ggplot2 is as a graphing package. After my experience I would use subprocessing into R rather than deal with rpy2 again

You're using an outdated version of R. You can use the R 4.0 from conda-forge channel, just create an environment with r-base and then install rpy2, the latest version compatible with R 4.0 is available through pip
conda create -n rpy2 -c conda-forge r-base
conda activate rpy2
pip install rpy2

Now you can use rpy2 with R 4.0 in your python code.
import rpy2.robjects as robjects
robjects.r('install.packages("BiocManager", '
           'repos="http://cran.r-project.org")')
robjects.r('BiocManager::install("seqLogo")')
robjects.r('library(seqLogo)')

Remember deactivate your env with conda deactivate.

Thank you so much this work! but when I do the same thing for "GenomicAlignments" package I get error:
raise embedded.RRuntimeError(_rinterface._geterrmessage())
rpy2.rinterface_lib.embedded.RRuntimeError: Error in (new("nonstandardGenericFunction", .Data = function (expr, envir = parent.frame(),  :
object 'GenomicAlignments' not found

I can't replicate you error. I did robjects.r('BiocManager::install("GenomicAlignments")') and then robjects.r('library(GenomicAlignments)') and it works.

I am unable to do this for pcalg (https://cran.r-project.org/web/packages/pcalg/pcalg.pdf). Can I see a minimal working example?


text:
How do I lift GWAS results to hg38?
There are some questions about lifting between reference builds, e.g. this one. But there doesn't appear to be a question about lifting a GWAS results file to a new reference build (except off-site). GWAS results files are not VCFs or plink files, so they aren't covered by the usual lifting tools (such as CrossMap or LiftOver). How do I lift a GWAS results file to hg38?
For instance, the Lee et al 2018 GWAS on education has this file:
user@computer:/home/user/Lee_2018$ head GWAS_EA.to10K.txt
MarkerName  CHR POS A1  A2  EAF Beta    SE  Pval
rs9859556   3   49455986    G   T   0.6905  -0.02901    0.00151 4.61e-82
rs7623659   3   49414791    T   C   0.3095  0.02899 0.00151 6.05e-82
rs11917431  3   49644012    C   T   0.6973  -0.0292 0.00152 7.40e-82
rs1873625   3   49666964    C   A   0.6973  -0.02916    0.00152 8.10e-82
rs11921590  3   49644193    T   C   0.6973  -0.02916    0.00152 1.17e-81
rs2352974   3   49890613    T   C   0.5 -0.02569    0.00141 6.35e-74
rs7029718   9   23358495    G   A   0.5646  -0.0244 0.00143 7.82e-65
rs7868984   9   23357826    C   T   0.4388  0.02437 0.00143 8.21e-65
rs11793831  9   23362311    T   G   0.4354  0.02418 0.00143 9.34e-64

Yes. LiftOver tools, if they work on VCF data (some only take BED as input) work on GWAS summary statistics. You could use
https://pypi.org/project/liftover/, for example.

Yeah, I figure out the simpler solution is to make a small VCF with an empty case, and then lift that to hg38 using liftover, then merge back the column of new positions.

I would advise against using dbSNP rsIDs to perform a liftover as they get updated over time. I recently wrote a software to liftover summary statistics and the documentation includes instructions that show you how to easily generate educational attainment summary statistics for GRCh38

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.

Once the plugins have been correctly installed, you could lift your GWAS results with a command like this:
bcftools +munge -f human_g1k_v37.fasta -C colheaders.tsv GWAS_EA.to10K.txt | bcftools +liftover -- -s human_g1k_v37.fasta -f Homo_sapiens_assembly38.fasta -c hg19ToHg38.over.chain.gz

It turned out to be easier than expected. You figure out which dbSNP version they used. dbSNP is the classification system for snps based on the rsid format. Pretty much every GWAS publishes their results with the rsids for the variants, as done above as well (MarkerName). So you download one of the dbSNP files from the official site. In this case, I am using version 142 since this was used by the 1000 genomes dataset.
After this, it is a matter of joining the tables. This can be done easily in R (you need about 15GB free memory), we load the dbSNP database into memory. Since we don't need all the columns, we skip some of them to save memory:
#library
library(tidyverse)

#cols
vcf_cols = cols_only(
  `#CHROM` = col_character(),
  `POS` = col_number(),
  `ID` = col_character(),
  `REF` = col_character(),
  `ALT` = col_character()
  )

#this is from the 1kg
#ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/other_mapping_resources/
dbsnp142 = readr::read_tsv(
  "ALL_20141222.dbSNP142_human_GRCh38.snps.vcf"),
  comment = "##",
  col_types = vcf_cols
  )

GWAS_EA = read_tsv(GWAS_EA.to10K.txt,
                   col_types = cols(chrpos37 = col_character(),
                                    chrpos38 = col_character()))

#make chrpos
dbsnp142 %<>% mutate(
  chrpos38 = str_c((CHROM %>% str_replace("chr", "")), ":", POS)
)

#count overlap
(GWAS_EA$MarkerName %in% dbsnp142$ID) %>% table2()

#mutate and merge
GWAS_EA = GWAS_EA %>%
  mutate(
    chrpos37 = str_c(CHR, ":", POS)
  ) %>%
  left_join(dbsnp142 %>% select(ID, chrpos38), by = c("MarkerName" = "ID"))

#update file on disk
GWAS_EA %>% write_tsv("GWAS_EA.to10K.txt")

The two new columns, chrpos37 and chrpos38, have the chr:pos variant names in both versions.


text:
Merging plink files with divergent allele codings
I have a number of plink files from different arrays (mixed Illumina and Affymetrix) that need to be merged. However, these use different allele codings, so trying to merge without fixing this causes an error (3+ alleles). One can fix this by using (on each file) --recode 12 to a PED file, and then converting back to BED (--make-bed) and then merging (--merge-list). Doing this, I was able to merge the files without explicit errors. However, this method is suboptimal, I think, because a minor allele in dataset 1 might not be the same one as in dataset 2, since this depends on the people in the dataset.
For instance, rs746969. Illumina 1MDuo:
10  rs746969    1.00698e+06 100008620   1   2

Affymetrix 6.0:
10  rs746969    0   100008620   A   G

Alternatively, I could keep the 12-encoded merged data, and attempt to convert back to ACGT. Which is the best approach?

Yes I would advise against doing the first option, since like you say, you don't know the minor allele is going to be the same between the different datasets. When merging two different datasets, it is very important that you make sure there aren't any strand flips and that the alleles line up. I would reccomend making sure that both datsets are in the ACGT format and then using something like genotype harmonizer to check for strand flips.


text:
Microarray Meta-Analysis and Cross-Platform
I have a microarray dataset with Agilent platform, all of the samples are cancerous, I need normal or control samples to compare with cancerous samples, can I combine control samples from another platform (for example Illumina) to Agilent platform? and then remove the batch effect (by ComBat function in the sva package).
I would be grateful if you could guide me to do it.

Hi there isn't sufficient detail really but I'll pop an answer below

In short the answer is no because you have to standardise between the arrays. Thus both arrays would need to run exactly the same control sample and you would then perform statistical standardisation (its a formal term), which is a varience measure where you look at the variation around the 'mean' for all samples in the array and here the 'mean' is the standard sample. If that can be done for both arrays at that point they can be combined.
The problem is intensity per run varies and you have to control against this

In this scenario you  can certainly combine samples by assessing whether the control samples between runs approximates to identity. What you could use is a uniform distribution to assess this, variance must MUCH smaller than the mean for the STANDARD sample(s). I don't really do uniform distribution statistics, so I'd need to research this, because 'uniform' means identity to standard deviations against the uniform distribution is a bit meaningless. However again the uniform distribution is denoted by a very variance in comparison to the mean. The statistic of dividing mean by standard deviation/variance from memory is called the Index of Association (I could be wrong). If this value is at an extreme you're good to go <<<1, if it is remotely close to 1 you're in trouble, if it is >1 you need to use standardisation.
Presumably you have several standards to ensure a mean is meaningful. A T-test and non-parametric equivalents is not stringent enough. If there is a difference assess a correction factor and rerun the uniform distribution check. The same correction will be applied to all samples.
The undisputed method to combine is standardisation, but accept in this scenario the differences between runs is potentially very small, so here what you are attempting is essentially to combine raw data or at least keep one of the runs as pure raw data.

Apologies for the delay in replying. Thank you for your valuable guide and useful information.
Can I combine control samples from the same platform (Agilent) to the desired dataset?

@MortezaHadizadeh I've answered above. If this answers your question please tick as 'accepted'.


text:
Converting a matrix array to matrix
I have a matrix in R. When I run class(my.matrix) it returns:
"matrix" "array"

Why do I get "array" when in fact I want a plain matrix? The same class output is returned when running as.matrix(my.matrix)).
Why is that?

From the return value of the class function, it looks like this already is a matrix. It just happens that it is considered by R to be both a matrix and an array.


text:
TCR-seq or scRNA-seq
I need idea, intuition, suggestion please
I have 10X scRNA-seq of PBMC (multiome's 3' poly-A capture), can I capture TCR/BCR from my scRNA-seq or I need TCRseq ? I googled and likely TCR and BCR immune profiling require 5' transcript capture and the V(D)J variable sequence is on the 5' end of the gene. So the 3' capture typically yields the C constant region of TCR and BCR transcripts. can you help me please? Thank you

The 3' gene expression protocol will capture TCR and BCR mRNAs but this may not be very helpful to you. As you already mentioned only the 3' end will be sequenced, which are the constant regions. With this data you may potentially be able to discriminate TRAC from TRBC1 and so forth.
However, in most cases researchers are interested in clonotypes (combinations of variable regions) and these you can only obtain by using the dedicated TCR or BCR kits from 10X. They contain specific sets of primers for these genes to cover the whole TCR or BCR gene body and not only the 3' end. This way you can see the different clonotypes in your population
The human TCR loci can be found for example at HUGO. The last gene at each locus is the constant and is identified by a C at the end. TRAC = T cell receptor alpha constant (I did not use the correct naming convention in the original answer, sorry)
B-cell receptors basically are membrane bound immunoglobins (IG). Therefore you will likely find the constant regions of the heavy chain in your data. The different heavy chains are what identifies the IG. For example IGHA1 is the heavy chain of an alpha IG. If the secreted and membrane bound form of the IGs is always the same per B-cell, unfortunately, is beyond my knowledge of B-cells.

Thank you I got a question, in this data base I could not find anything corresponding to TCRVδ2 , you think in which part I can find this gene?

The nomenclature used by immunologists may not match the official gene names (which are in HUGO). The TCR gamma variable genes are named TRGV followed by a number. Problems with different names for the same gene / protein are all over the place and you will sometimes have to search several data bases (HUGO, UniProt, Entrez ... ) or the literature for an answer

According to our experience 3’ sequencing does not produce enough data to recover TCR and BCR repertoire. But we have some experience recovering repertoires from non enriched 5’ 10x RNA-Seq GEX data.
With MiXCR we have specific presets that allow analyse 10x data in one single line command:
mixcr analyze 10x-5gex-full-length \
 --species hsa \
 sample_R1.fastq.gz \
 sample_R2.fastq.gz \
 sample_result

Also feel free to check out our documentation portal:
https://docs.milaboratories.com/mixcr/reference/overview-built-in-presets/#10x


text:
Help with Seurat QC ambiguity
I have four PBMC samples from 10X scRNA-seq
> cancer
An object of class Seurat
36601 features across 18338 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

And this is Seurat QC plot

If I am not wrong, samples were of relatively low quality, with gene expression data revealing the presence of mitochondrial genes, as well as MALAT-1, which are suggestive of poor sample quality (dead/dying cells).
Any way by this plot I filtered cells to remove data that had more than 20 "percent" mitochondrial expression , > 2000 features or < 100 features
cancer <- subset(cancer, subset = nFeature_RNA > 100 & nFeature_RNA < 2000 & percent.mt < 20)

By this I lost most of cells
> cancer
    An object of class Seurat
    36601 features across 6,883 samples within 1 assay
    Active assay: RNA (36601 features, 0 variable features)
    >

Please, can somebody experienced in scRNA-seq inspect this plot and my thresholds and tell me if I am wrong, any suggestion, if losing this number of cells is normal
Thank you

Answer from @haci, converted from comment:
I think you should not be using a cut-off higher than 20 or 25 for the reasons I explained above.
I would advise against using cells with more than 20-25% percent.mito. This metric is used as a proxy for cell death and there is not much point keeping these cells if the cell death is so severe to the point where let's say half the reads map to the mitochondrial genome (with a 50% cut-off). There will be many caveats one being apoptotic cells from different lineages clustering together, not interesting unless you are studying apoptosis related phenomena. In my opinion 6000 good quality cells are better than 10000 apoptotic cells.

From the graphs, you're losing most of the cells through the percent.mitochondrial threshold. The distribution is so broad on that mitochondrial plot that I probably wouldn't filter for mitochondrial expression (except possibly at ~95%). If I did have a lower threshold, I'd put it at around 50% rather than 20%, because that's where there's a change in density on the mitochondrial plot.
This suggests that there's a very large mitochondria-rich population in the PBMC dataset. I don't have a good understanding of biology and cell populations, but if I were presented with these data I'd probably try clipping out those cells with 65-95% mtDNA and see if there is a strong enrichment for markers that correspond to any particular cell population (e.g. platelets).


text:
Filtering cells from values in metadata
I metadata slot of Seurat object I have mapping score of each cell to a reference PBMC data like
pbmc <- subset(object, subset = object@meta.data$mapping.score > 0.5)
Error in FetchData(object = object, vars = unique(x = expr.char[vars.use]),  :
  None of the requested variables were found:
>
>
>

I want to only keep cells with mapping score > 0.5
Any help please?

your_new_object <- subset(your_object, subset = mapping.score > 0.5)


text:
which marker identifications method is recommended when single-nuclei RNA seq datasets are integrated across regions
I'm analyzing Sn-RNA seq datasets which include different brain regions and spinal cord in patients and normal controls. I would have different comparison within and between individuals. I've got a question which I could not find the answer. First, I know that in the clustering step FindConservedMarkers is the best choice when you are comparing controls & patients. I was wondering does that work the same when I want to compare two different regions like the brain and spinal cord? if not what do you recommend for this comparison? as there might be different cell population in two different regions and those regions are not the same

Hi @Paria I can't quite see how this is a question, its either a 'yes'  or a 'no;' (I'd steer towards yes). Maybe if you explain how you arrived at the conclusion that FindConservedMarkers was 'best choice' an informed point of view can be forwarded.

Hi @Michael. I am saying based on the Satija lab workflow. I have single-nuclei 10X genomics datasets. FindConservedMarkers code in Seurat helps to identify the shared gene expression across condition and cell-type identification. In my case I'm comparing two different region within individuls (same condition different regions) and I am asking wether this work the same in my case?


text:
Compute the significance of the overlap between 2 or more gene sets
I was able to compute the significance of the overlap between 2 gene sets using the cdf function of scipy hypergeometric distribution.
I wish to be able to perform the same calculation for more than 2 gene sets; should I use the multivariate hypergeometric distribution cdf function for that?
Are there any websites that provides the same calculations over gene sets so I can validate my results?

Not so trivial, if you want to use a distribution. https://stats.stackexchange.com/questions/52004/overlapping-sets-3-significance-test

Easiest way is to simulate the data, and see how many times you end up with an intersection set larger than what you see

This is a shot at it, first an example dataset:
import matplotlib.pyplot as plt
import numpy as np
import functools
from matplotlib_venn import venn3

# define universe
uni = ["gene"+str(i) for i in range(1000)]
# some overlap
gs1 = uni[250:300] + uni[900:950]
gs2 = uni[:300]
gs3 = uni[250:500]

The ever amazing venn diagram:
venn3([set(gs1),set(gs2),set(gs3)],set_labels=["gs1","gs2","gs3"])

Then a function to draw a set with length equivalent of each set, randomly from the universe and find length of intersection (all 3):
def sim_intersect(uni,set_lengths):
    randomsets = [np.random.choice(uni,n) for n in set_lengths]
    return len(functools.reduce(np.intersect1d,randomsets))

We run this 1000 times:
permuted_values = [sim_intersect(uni,[len(gs1),len(gs2),len(gs3)]) for i in range(1000)]

plt.hist(permuted_values,bins=range(50))

The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this post:
(sum(np.array(permuted_values)>obs_n)+1)/(1000+1)
0.000999000999000999

How would a population size parameter fit in this approach? Doesn't the p value should be affected by the total 'available genes' ?
For example,. if I have two gene sets A, B.
|A| = 250, |B| = 220, |A&B| = 65 As far as I understand,
the p value should change dramatically if the population size is 300 or 30,000.
Could you please explain where/how does it fit in this approach?

Yes you are right, the total available genes would matter. In the code above, you would define the set of available genes as uni . It is not like a hypergeometric where you input a number of genes parameter. You simply simulate the results and see how likely you end up with a result as extreme as yours

easiest thing to do is to re-run the above simulation with a small "universe"


text:
How to refine this distribution comparison script
My data look like this

They are compound similarity estimations(the whole file is around 10GB). What I am trying to achieve is to compare the similarity distributions of each compound using the Kolmogorov-Smirnoff test so as to use the data for AIC calculation.
This is my code:
import pandas as pd
from scipy.stats import ks_2samp
import re

with open('ecfp4_file.csv', 'r') as f, open('Metrics.tsv', 'a') as f_out:
    f_out.write('compound_1' + '\t' + 'compound_2' + '\t' + 'Similarity' + '\t' + 'KS Distance' + '\n')
    df = pd.read_csv(f, delimiter = ',', lineterminator = '\n', header = 0)

    for m in range(10, 99):
        dc1 = []
        dc2 = []
        m1 = m/100
        m2 = m/100+0.009
        df2 = df.loc[(df['ECFP4'] >= m1) & (df['ECFP4'] < m2)]          #Work on different thresholds each time

        for i in range(0, df2.shape[0]):
            if len(df2.index[df2['id1'] == df2.iloc[i, 0]].tolist()) >= 5:
                dc1.append(df2.iloc[i, 0])
                dc2.append(df2.iloc[i, 1])
        for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            f_out.write(str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

which so far has taken 10 hours to run and it has not yet passed the 0.1 threshold, so it is highly unpractical. Does anyone know of a better way to do this? Or can help me refine my code for faster and better results?
I also tried to refine it and run it in R, but I get a vector memory exhausted (limit reached) error when I run it
install.packages("ggplot2")
library(ggplot2)
PS22_ECFP4 <- read.csv("./PS22_ECFP4.csv")

df <- data.frame(matrix(ncol = 4, nrow = 0))
cols <- c("id1", "id2", "Similarity", "KS Distance")
colnames(df) <- cols

for (i in PS22_ECFP4) {
  sim_df <- PS22_ECFP4[PS22_ECFP4$id1 == i,]
  dist1 <- sim_df[ , 3]
  for (k in 1:nrow(sim_df)) {
    sim <- sim_df[k , 3]
    sim_r <- floor(sim * 100) / 100
    comp2 <- sim_df[k, 2]
    comp2_df <- PS22_ECFP4[PS22_ECFP4$id1==comp2,]
    dist2 <- comp2_df[ , 3]
    if (length(dist1) > 1 & length(dist2) > 1) {
      ks <- ks.test(dist1, dist2)
      df[nrow(df) + 1,] = c(i, comp2, sim_r, ks$statistic)
    }
  }
}

write.csv(df, "./Metrics.csv", row.names = FALSE)

I am not sure how useful this answer is. But may be you could try this.
For your python script:
You use f_out.write within a loop. I think usually this I/O operation is resource intensive. I guess the speed would improve significantly if you store the string in a variable and write the output only once outside the for loop like this:
out_str = ''
for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            out_str +=  (str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

f_out.write(out_str)

Regarding the memory error in your R script, are you using 32bit build of R?
Please see this:
http://web.mit.edu/~r/current/arch/i386_linux26/lib/R/library/base/html/Memory-limits.html#
They say "Under most 64-bit versions of Windows the limit for a 32-bit build of R is 4Gb: for the oldest ones it is 2Gb. The limit for a 64-bit build of R (imposed by the OS) is 8Tb."
Conclusion: Please try using your python script, but use the write operation outside the loop.

If I was doing something like this in R (10G file), I'd use data.table
I can't quite work out what you're R code is doing (does it do what you want it to?)
I think the below might do what you want
library(data.table)
dat <- fread("./PS22_ECFP4.csv")
dat[,Similarity:=floor(ECFP4*100)/100]
dat[,'KS Distance':=lapply(1:nrow(dat),function(i)
  ks.test(dat[id1==dat$id1[i],ECFP4],dat[id2==dat$id2[i],ECFP4])$statistic
)]
dat[,ECFP4:=NULL]
fwrite(dat,"./Metrics.csv")

Edit:
The script loops through the data, and for each row creates two subsets of the data matching id1 and id2. Then does the ks.test on the similarity columns of the two subsets.


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

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

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

in genomics, organisms have a single reference genome which is updated( or not) over a period of time. E.coli has different strains and this means there might be reference genomes for each of the strain. If there is no reference genome available the its likely that in constructing the sequence of those strains they might have used any of those two as their reference. This information you can find by reading the paper that published those strains. Therefore I will advice you pick any of those two available reference genome as your reference and work with and make sure you indicate which one you used in your report or document.


text:
Standard Way to Preprocess Gene Expression?
I am trying to collect gene expression data for the point of fitting gene regulatory networks. My background is primarily in computer science and so I am finding the biological literature a bit difficult to penetrate.
Specifically, I want to fit two regulatory networks for several cancers in TCGA from tumor and normal samples. I have queries using the R package TCGABiolinks and was able to follow the tutorial here up to the normalization part.
So I have three questions;

Is there a standard way to normalize gene expression counts? If so, what is it, and is it available as a library in R, python or some other language?
If not, can I just center and scale (z-transform) the samples? I am not interested in differential expression, only quantifying expression from counts.

Question is, what do you mean by 'quantifying expression from counts'? TCGA has a standard normalization available, which is FPKM. It normalizes expression by fragment-size and library depth. Some people prefer TPM to FPKM though (it alters the order of operations), so you may want to download the count data and do a TPM normalization using edgeR.  For visualization, for TCGA, I personally tend to log2(FPKM) and row-center for heatmaps.

okay so is the FPKM normalization what is given in the scaled_estimate design matrix?

Not sure what that matrix is... in TCGAbiolinks or any other way to download from the GDC, you'll have FPKM. It is not scaled though, so probably no.

First off - do not use FPKM for across-sample normalization. This is such a common issue that I have a prepared text for it:

Why not to use FPKM

They have a poor sensitivity when looking at low expressed genes [1, 2]. Short transcripts/genes may have falsely over-estimated expression values [5].
They may not normalise adequately across samples [3]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in [4]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself [6].
They have a high False Positive Rate in differential expression analysis [1, 2]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.

On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected [7]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.

References
[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package:
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839

What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function varianceStabilizingTransformation() from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions voom, which is another good alternative, as is limma. All of these are available on Bioconductor.


text:
Basic Docking Work Flow
I'm a computer scientist who is now getting into designing new molecules, with the goal that they dock well to a given target. While I'm used to coding up machine learning algorithms in python, I'm finding things a bit murky in terms of the how-to of docking. I lack pretty much any (medicinal) chemistry/biomedical engineering background and so understanding how to get things to work has been a real struggle.
In particular I've been looking at two githubs: this one using smina and this one using glide. Recently, my collaborators have steered me toward rDock/RxDoc.  Near as I can tell the basic ingredients for ANY docking software are:

a target file (rna or protein) in a .pdbdt,.pdb, or .mol2 file with
it docked to any ligand
a ligand file (mol2)
a target file with the ligand removed

But where I get lost is in any intuition on how to prepare these files, how to find them, and any intuition at all about how to make these work. I've seen suggestions of openbabel, etc, but there seem to be no tutorials online.
So, my question is: Suppose I have a molecule in rdkit and a protein or an RNA that I want to dock it to, accessed from the protein databank, say. What is the ``right'' way to do the following:

Download the receptor from pdb
Process the receptor to get file of the receptor with and without ligand
Process a molecule from rdkit so that the docking software can dock a (new) molecule to the receptor and give a docking score?

Ideally this should work for `any' receptor, but if the steps above need to be based on the specifics of the receptor, that would be helpful for me to know.
I'm shocked there isn't a beginners guide on this out there that's more easy to find, and I worry about going down this route without any intuition.

The basics of docking are explaining in several textbooks (eg. Essentials of computational chemistry). The science is rather important as there are lots of complications as protein are not solid blocks, water can form tight interactions (cf implicit water), Lipinski rules, the candidate compound needs to block the protein s activity (cf box size), Gibbs free energy (“score”) isn’t the sole value to sort on etc etc.

Regarding what to use, the Capri competition is a great way to see how programs compare. Glide scores well but it’s not open source. Autodock4 and Autodock Vina are very popular and highly documented —they are different in approach and the former is more adaptable and the latter is (misleadingly) simple. The papers describing these give a lot of info…

not sure I remembder it right, you put your target in a box , add the ligand somewhere in the box , minimize everything to find the most probable binding site ?, or you force the ligand in one particular site of the target ?