How to map PDB chains to Uniprot IDs using API services

I found the answer and I thought it may be useful for others. we can use pypdb package for that purpose (it works in unix systems):

import pypdb

all_info = pypdb.get_all_info('1kf6')
print(all_info)

the output is as follows:

{'polymer': [{'@entityNr': '1', '@length': '602', '@type': 'protein', '@weight': '66057.6', 'chain': [{'@id': 'A'}, {'@id': 'M'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase flavoprotein subunit', 'accession': {'@id': 'P00363'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE FLAVOPROTEIN'}, 'enzClass': {'@ec': '1.3.5.4'}}, {'@entityNr': '2', '@length': '243', '@type': 'protein', '@weight': '27021.9', 'chain': [{'@id': 'B'}, {'@id': 'N'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase iron-sulfur subunit', 'accession': {'@id': 'P0AC47'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE IRON-SULFUR PROTEIN'}, 'enzClass': {'@ec': '1.3.5.1'}}, {'@entityNr': '3', '@length': '130', '@type': 'protein', '@weight': '14898.8', 'chain': [{'@id': 'C'}, {'@id': 'O'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit C', 'accession': {'@id': 'P0A8Q0'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 15 KDA HYDROPHOBIC PROTEIN'}}, {'@entityNr': '4', '@length': '119', '@type': 'protein', '@weight': '13118.9', 'chain': [{'@id': 'D'}, {'@id': 'P'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit D', 'accession': {'@id': 'P0A8Q3'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 13 KDA HYDROPHOBIC PROTEIN'}}], 'id': '1KF6'}

This output gives all the information related to all chains and we can simply find uniprotkb AC of each chain in this output. I also found another answer for this question. This answer can be found here and the following code:

import requests
from xml.etree.ElementTree import fromstring

pdb_id = '1kf6.A'
pdb_mapping_url = 'http://www.rcsb.org/pdb/rest/das/pdb_uniprot_mapping/alignment'
uniprot_url = 'http://www.uniprot.org/uniprot/{}.xml'

def get_uniprot_accession_id(response_xml):
    root = fromstring(response_xml)
    return next(
        el for el in root.getchildren()[0].getchildren()
        if el.attrib['dbSource'] == 'UniProt'
    ).attrib['dbAccessionId']

def get_uniprot_protein_name(uniport_id):
    uinprot_response = requests.get(
        uniprot_url.format(uniport_id)
    ).text
    return fromstring(uinprot_response).find('.//{http://uniprot.org/uniprot}recommendedName/{http://uniprot.org/uniprot}fullName'
    ).text

def map_pdb_to_uniprot(pdb_id):
    pdb_mapping_response = requests.get(
        pdb_mapping_url, params={'query': pdb_id}
    ).text
    uniprot_id = get_uniprot_accession_id(pdb_mapping_response)
    uniprot_name = get_uniprot_protein_name(uniprot_id)
    return {
        'pdb_id': pdb_id,
        'uniprot_id': uniprot_id,
        'uniprot_name': uniprot_name
    }

print (map_pdb_to_uniprot(pdb_id))

The output of this code is as follows:

{'pdb_id': '1kf6.A', 'uniprot_id': 'P00363', 'uniprot_name': 'Fumarate reductase flavoprotein subunit'}


Difference between samtools mark duplicates and samtools remove duplicates?

samtools rmdup and samtools markdup -r do the same thing. Without the -r flag samtools markdup only flags the duplicates.

You'll have to run samtools fixmate -m and sort the output to add ms and MC tags to prepare the file for markdup..

According to the current documentation, rmdup is obsolete. Please do not use rmdup.


Run kallisto iteratively across many samples

Running the jobs in parallel

Iteratively is one solution, but since you're not using bootstrapping, which can be run multithread (and is not needed if you're taking the popular tximport route) you can do better using GNU parallel, this will enable you to run as many jobs as you have execution threads on your system simultaneously.  Note that Kalliso only runs multithreaded with -t if you're using bootstrapping.

You'll need to make the index first:

kallisto index --make-unique -I index.fa.idx index.fa

Next lets define the input FASTQ files and the sample IDs you want to use using a tab-delimited input file like so:

sample_1    sample_1_R1.fastq.gz
sample_2    sample_2_R1.fastq.gz
sample_3    sample_3_R1.fastq.gz

Here the first column is the sample ID which will be used for output and the second separated by a tab character is the input file.

You can produce this using a text editor but as a short-cut you can feed the list of fastq.gz you have into a list via:

ls -1 *.fastq.gz > editme

Next we create a list of jobs to run on our system using a utility script Make_job_list_kallisto.sh which consumes our tab delimited input file list above.  I've used your -l and -s parameters below, note that with paired end data normally these are estimated for you.

#!/bin/bash -eu

[ $# -ne 2 ] && { echo -en "\n*** This script generates jobs for GNU parallel. *** \n\n Error Nothing to do, usage: < input tab delimited list > < output run list file >\n\n" ; exit 1; }
set -o pipefail

# Get command-line args
INPUT_LIST=$1
OUTPUT=$2

# Set counter
COUNT=1
END=$(wc -l $INPUT_LIST | awk '{print $1}')

echo " "
echo " * Input file is: $INPUT_LIST"
echo " * Number of runs: $END"
echo " * Output job list for GNU parallel saved to: $OUTPUT"
echo " "

# Main bit of command-line for job
CMD="kallisto quant -i index.fa.idx --single -l 200 -s 0.1"

# Main Loop
[ -e $OUTPUT ] && rm $OUTPUT
while [ $COUNT -le $END ];
do
    LINE=( $(awk "NR==$COUNT" $INPUT_LIST) )
    # Make file list
    echo "Working on $COUNT of $END Sample ID: ${LINE[0]}, Files ${LINE[@]:1}"
    echo "$CMD -o ${LINE[0]} ${LINE[@]:1}" >> $OUTPUT
    ((COUNT++))
done

We use the script like so:

 Make_job_list_kallisto.sh input_file_list.txt job_list.txt

Now we have a job list for GNU parallel we can run our jobs in parallel like so:

parallel --progress --jobs 4 --joblog kallisto_joblog.txt < job_list.txt

This will run 4 of the kallisto jobs simultaneously (Increase this if you have more threads - don't mind consuming all the IO/CPU threads on your system) producing a nice progress report of those completed, additionally a job log will be written to the file kallisto_joblog.txt.

For paired end data I already have these scripts + similar instructions on GitHub which people might find useful.


how to solve “dim(X) must have a positive length” at running ComBat function in R

batch = Pheno_LMS$batchId should be batch = Pheno_LMS$BatchId, if you fix that then at least the code you showed works just fine.


Bit of an old post sorry, but I think you might be getting a little confused here. The nf-core pipelines come with their own software already, so you don't need to install or bootstrap optitype yourself. If it's used in the pipeline, it comes in the pipeline container(s) (or conda environment or whatever packaging system you use). So you should never have to touch docker://fred2/optitype.
When you run the pipeline with the -profile singularity argument, it will automatically pull the image from Docker Hub and convert to Singularity for you (some of the newer pipelines pull pre-build Singularity images directly). You can also use the nf-core download hlatyping --singularity command from the helper tools to pull the image ahead of time.

Docker images for this pipeline: https://hub.docker.com/r/nfcore/hlatyping
Documentation for this pipeline: https://nf-co.re/hlatyping/
nf-core/tools download docs: https://nf-co.re/tools/#downloading-pipelines-for-offline-use

I hope this helps! Come and join the nf-core Slack organisation if you need more help. We have a channel for each pipeline, including #hlatyping.


Error while performing "lift-over" of genomic coordinates between two species using Ensembl Compara API

While I wouldn't consider this as the best answer or solution - I decided to use Ensembl REST api to extract the coordinates of alignment blocks.

The quick fix - in python looks like this:

#This is the main function which uses Ensembl REST API to access the Enesembl Compara database
#And retrive the alignment blocks which correspond to the MM10 standard.
#The function accepts chromosome, start and array and returns the alignment blocks
#reported by the API. Since the blocks can be on different chromosomes/and/or/contigs, it
#returns the whole output in form of dataframe - containing chromosome, start and stop of the
#alignment blocks corresponding to input and all the other alignment blocks corresponding the
#16 mice species, if available.
#!!!!!!!----------Caution start----------!!!!!!!
#The API is a bit buggy, so it is advised not to request for specific species.
#!!!!!!!-----------Caution end-----------!!!!!!!
def liftover(chrom, start, end):
    mmStart, mmEnd, mmchr = [], [], []
    caStart, caEnd, cachr = [], [], []
    #ryStart, ryEnd, rychr = [], [], []
    server = "https://rest.ensembl.org"
    ext = "/alignment/region/mus_musculus/" + str(chrom) + ":" + str(start) + "-" + str(end) + ":1?species_set_group=murinae;method=CACTUS_HAL;content-type=application/json"
    headers={ "Content-Type" : "application/json"}
    r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
    if(r.ok):
        json_string = json.dumps(r.json())
        jData = json.loads(r.content)
        results = r.json()
        print("The response contains {0} properties".format(len(jData)))
        print("\n")
    else:
        r.raise_for_status()

    for element in jData:
        element.pop('tree', None)

    for element in jData:
        for key, value in element.items():
            for z in value:
                if (z["species"]=="mus_musculus"):
                    mmStart.append(z["start"])
                    mmEnd.append(z["end"])
                    mmchr.append(z["seq_region"])
                if (z["species"]=="mus_musculus_casteij"):
                    caStart.append(z["start"])
                    caEnd.append(z["end"])
                    cachr.append(z["seq_region"])
    if (len(mmStart) != len(caStart)):
        print("Length mismatch between species blocks - proceed with caution!")
    return [min(mmStart), max(mmEnd), min(caStart), max(caEnd), len(mmStart)]

Again, it depends on how alignment blocks are reported - and requires further digging to resolve which mm10 alignment block resembles CAST alignment block.


Tool to show DNA sequence and allowing upload of own graph data

You can use UCSC genome browser, which supports custom genome sequence and annotation. Specifically, one needs to create a Track Hub to hold genomes and annotations.

Here is an example:

  http://genome.ucsc.edu/cgi-bin/hgTracks?hubUrl=http://zlab-trackhub.umassmed.edu/yfu/insects/hub.txt&db=hub_129615_cl_v0.8.2

And here is the corresponding hub.txt file holding the metadata, and URLs to actual fasta, bigWig and bigBed files.


What's the most efficient file format for the storage of DNA sequences?

Genomes are commonly stored as either fasta files (.fa) or twoBit (.2bit) files. Fasta files store the entire sequence as text and are thus not particularly compressed.

twoBit files store each nucleotide in two bits and contain additional metadata that indicates where there's regions containing N (unknown) bases.

For more information, see the documentation on the twoBit format at the UCSC genome browser.

You can convert between twoBit and fasta format using the faToTwoBit and twoBitToFa utilities.

For the human genome, you can download it in either fasta or twoBit format here: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/


How can I compute gene expression for a set of RNA reads?

It's unclear if your paired-end reads are actually in fasta format, I'll presume that they're in fastq instead.

The easiest tool to use is salmon, which nicely deals with things like multimapping. If you're trying to judge the quality of an assembly, then I recommend having a look at transrate, which uses some related methodology for assessing contig quality.


Voom transformation of RNA seq raw counts data

Try to use DElist() function before you transform, and also make rownames first.

counts <- prac_count_10[,-1]

rownames(counts) <- prac_count_10[,1]

DGE1 <- DGEList(counts)

And then continue with your voom...


Profile of conservation between bacterial sequences and human protein

I'm not sure what you mean by it being a mess? That looks like a pretty good alignment to me, and it most certainly has aligned all the sequences. You are nearly always going to have gaps (see my MSA below, which are all paralogs from the same genome so they couldn't really be more closely related, and yet, there are gaps). That comes with its own set of problems mind you, as you need to decide how you're going to deal with them.

You also don't necessarily have to do a profile alignment. Those sequences don't look massively divergent to me, so straight forward sequence alignment would probably give you more or less the same result.

  My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

This doesn't really make sense. You should really be searching for the number of BLAST hits with a particular domain, if you want to infer homology.

If you want a per-position visualisation of the conservation, you could look at the shannon entropy for each column and plot that. I wrote a script to do just that a little while ago: https://github.com/jrjhealey/bioinfo-tools/blob/master/Shannon.py

Just beware it's not super well tested yet. Feed an MSA in with as many sequences as you want to analyse, but you'll have to have identified the sequences and done the alignment first.

For example, given this MSA:

    16    149
PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTG
PLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTG
PAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTG
PLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTG
PAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTG
PLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTG
PLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIG
PAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVG
PAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVG
PAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PLT_02568  MSTTVDQIAV QYPIPTYRFV VTVGDEQMSF QSVSGLDISY DTIEYRDGIG

           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINITLRKGV FSGDTKLFDW LNSIQLNQVE KKDISISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NHYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NWFKMPGQSQ LVNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NWFKMPGQSQ STNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPSITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD

           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL ISPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TAPAFDATSN EVAVQEISLT ADRVTIQAA
           ETGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EVGTEILMTW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EAGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVMIQAA
           ETGSQILMTW NVANAFPTSF TSPSFDAASN DIAIQEIALV ADRVTIQAP
           EAGTEILMTW NVANAFPTSF TSPSFDATSN EIAVQEIALT ADRVTIQAA
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLT ADRVIMQAV
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLM ADRVIMQAV
           DTGSEVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           DTGSKVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEMSLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH

You'd get this plot:

  I would like to have a graph with my protein on the x-axis and a plot
  like the regions of high conservation for the multiple alignments,
  except with the spikes corresponding to a high number of BLAST hits.
  This would allow me to identify regions of my protein that have higher
  sequence similarity to bacteria than others.

I'm not sure your logic is quite right here though. Blast won't give you hits depending on a particular position. It's a local aligner, so it'll just return you hits where at least some part of your query matches at least some part of another.

What you could do is take the logic in the script above, and just use a different metric. For example, perhaps you could count the proportion of sequences which have the most common amino acid at a given position within your MSA. That would be fairly crude though.

As you say in the comments,

  My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

your original alignment will show you this intrinsically, if only you include the sequences of all the BLAST hits in the first place. Thus your work flow will be:

Blast sequence of interest.
Download all/as many hits as you want (bear in mind the E-value/bitscore and the number of hits you get. It might only be a few dozen, in which case you can use the lot, but if not, just take all the hits below a certain cut-off.)
Align all the sequences.
Look at the column scores for the whole MSA. You can use whatever metric of conservation you like really. Might be as simple as proportion of sequences with the most common residue, or something more complex like Shannon entropy (though as you can see in the graph above Shannon entropy can be kinda noisy) etc.


Finding drugs based on gene expression data

It would be great if you can provide more information about the gene expression data set. If you have response to single drug, the general approach is to find the deferentially expressed genes with tools like DESeq2, or make a heat map based on the expression of your listed genes to see the drug effect before and after the treatment. But, if you want to study the effects of combination of drugs it is necessary to consider the additive and synergistic effect of drug combinations.


Reordering scaffolds according to a reference without a genetic map

Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need.

Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a reference and/or b) rearrange contigs to improve assembly quality.

I use the tutorial as above which uses the GUI for Mauve. I


Recommended workflow for RNA long read Kallisto-like TPM estimation?

Stringtie now does long reads too. We've found flair to be quite difficult to install and to crash often. That said, stringtie isn't very good about showing you its working (e.g. I have no idea how to get it to give the assembly it uses to figure out all the isoforms). Its values for TPM/FPKM also seem a little weird, in that they're perfectly correlated with each other, and for transcripts they're perfectly correlated with coverage too.
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1910-1
https://github.com/gpertea/stringtie
TALON is another tool that does isoform assembly. I'm not sure whether it does quantification though.
https://github.com/mortazavilab/TALON#


Filter out outliers of the scRNA-seq (heterogenous cells)

From what I known, there is no clear consensus in the field and it depends on the type of cells you are interrogating.

However, if mitoRNA/endogenousRNA ratio does not fit for your purposes, other option is to check at the total number of genes/transcripts detected in each cell. In this way, you can filter out cells that have considerable less genes/transcripts detected than the rest since this may be indicative of bad quality cells for any reason (apoptosis, RNA degradation, sequencing itself, etc). For instance, to put a threshold of >1000 genes/transcripts in a cell to be considered for further analysis.

Consider also to take a look at the number of mapped reads for each cell, as cells with low mapped reads respect to the rest could be potentially problematic.


assigning different characters to present and absent genes from blast output using perl or python

I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk:

$ awk 'NR==FNR{a[$1]=0; next}
       {
        if($2 in a){
            a[$2]=1
        }
       }END{
        for(i in a){
            print i,a[i]
        }
       }' pangenome.txt blast.individual1.txt
Ha1_00044161 0
Ha4_00025807 1
Ha11_00031248 0
Ha16_00021664 1
Ha8_00042762 1

The above is just a one-liner, I expanded it for clarity. You can run it on as many files as you like:

awk 'NR==FNR{a[$1]=0; next}{if($2 in a){a[$2]=1}}END{for(i in a){print i,a[i]}}' pangenome.txt blast.individual*.txt

For many large files, the fastest approach would probably be to just grep the names (unlike the awk, this will also match if the target gene is in the 1st field, but that's probably not an issue for you):

$ grep -wof pangenome.txt blast.individual1.txt.orig
Ha16_00021664
Ha8_00042762
Ha4_00025807

That will give you the list of genes found in the target blast files. To get the 1 and 0, you can do:

$ sort <(grep -owf pangenome.txt blast.individual1.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'
Ha1_00044161 0
Ha11_00031248 0
Ha16_00021664 1
Ha4_00025807 1
Ha8_00042762 1

For multiple files, just run it on all of them:

$ sort <(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'


The Dockerfile is generated by galaxy-utils, so I'm not entirely sure what it looks like. However you should be able to generate such a "mulled" container as documented here. In short, something of the form:
mulled-build build 'bioconductor-csaw=someversion--somebuildnumber' 'someotherpackage=someversion--somebuildnumber' ...

That probably requires setting the CONDA_IMAGE environment variable. Bioconda currently sets that to quay.io/dpryan79/mulled_container:latest.


How to score how densely bases are distributed along a given oligo?

You could count the number of transitions between an A and a 'not A', and divide this number by the total number of As.  You might also need to scale by the sequence length, if you have different length sequences that you want to compare.

Edit: also, there might be a correction to make if one or both ends of the sequence are an A


How to filter out partially similar strings from two lists and make one list?

You really need to learn some python :)
with open('/path/list2.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip())

with open('/path/list1.txt', 'r') as fin:
    with open('/path/list3.txt', 'w') as fout:
        for line in fin:
            if line.strip().split(' ')[0] in to_keep:
                fout.write(line)


Can I stop my nanopore sequencing run if there are no more reads being produced?

Yes, you can click "stop acquisition" and your run won't be negatively affected. All of the reads are saved as they are generated. I am not sure how this will impact live basecalling though if that is something that you do.


Since this question's been asked there's been a (what looks like a first useful?) release of Drug Ontology DRON. Perhaps worth checking again?


What does "ALT+ cancer cells", the plus sign mean?

Here are a couple of papers (A and B) that explain ALT.

  Telomere maintenance is a hallmark of cancer as it provides cancer cells with cellular immortality. A significant fraction of tumors uses the alternative lengthening of telomeres (ALT) pathway to elongate their telomeres and to gain an unlimited proliferation potential. Since the ALT pathway is unique to cancer cells, it represents a potentially valuable, currently unexploited target for anti-cancer therapies.
  B

The ALT+ denotes cell lines that are positive for the ALT pathway.


Merging transcriptomes coming from different experiments

As mentioned in the comments, your biological question is going to be vitally important in determining the best way to tackle this issue. However, there are some things that you should generally keep in mind.

Most de novo assemblers will perform kmer normalization. This means that the merged RNA-Seq data will be normalized across all libraries. This is almost certainly not what you want because you'll be randomly selecting reads from all libraries without regard to their treatments. So, if you're going to allow for kmer normalization, you should normalize prior to merging the libraries so that you have equal representation of the reads by library. If you don't normalize, you're going to have a very fragmented assembly, and it will take a long time to generate.
De novo assemblers using short reads assume all identified exons with junction support may be used to generate splice variants with equal probability. I'm not entirely sure how to best phrase that, but here's an example. Suppose in a cold-treated library that Gene A generates a transcript using exons 1-2-3, but a heat-treated library generates transcripts using 3-4-5 and 2-3-5. It may be that merging the libraries results in transcripts generated containing various permutations of those exons (e.g. 1-2-3-5, 2-3-4, etc.) that aren't generated in any of the biological settings but can no longer be distinguished. This will lead to false positive transcripts which will detrimentally affect expression quantification.
Merging assemblies is difficult. So, to be true to each treatment, it makes more sense to generate individual assemblies (which will also result in library-specific kmer normalization). StringTie essentially does this using reference genomes. You assemble using each library (all bio-reps for a treatment), and then merge after the fact and re-quantify. This is the best approach. However, it can be difficult and computationally taxing. Some programs like CD-HIT (specifically cd-hit-est) are typically used in these cases. Keep in mind that CD-HIT will drop sequences with high similarity, not merge them. This may lead to loss of some uniquely assembled sequence. Other approaches involve using programs like CAP3 (a DNA assembler) to treat the transcripts like longer reads for assembly. Depending on the parameter choice, this might work, but you will risk generating chimeric transcripts.

How you move forward with the assembly process will again depend on your biological questions. I would tend to promote that you tackle the questions one at a time rather than using a single transcriptome to address them all, but there may be circumstances in which merging assemblies makes more sense.

In your specific case, I would perform multiple assemblies using a variety of the approaches described above and then annotate the assemblies (ex. using Trinotate) to see how their quality metrics compare. I would discourage generating a single transcriptome and then re-assembling individual transcriptomes from that assembly -- the single transcriptome may have misassemblies that would then propagate to the individual assemblies, and this process would generally be messy/difficult.


If need to get msa between every read and ref in the specific region ,
you can try this https://github.com/orangeSi/bam2msa.git

bam2sam ref.fa input.bam --regions chr1:200-300


Is there a way to get the code from "github.mit.edu"?

It is very strange that they put their code somewhere cannot be easily accessed. If they just put it under a normal Github repo, there would not be such a problem. One way to get the code is simply to email the corresponding author. People are usually responsive if you are interested in their work.


PCA plot did not work in single cell RNA-seq

The error is pretty obvious. There is no PCA data in the respective slot of the SingleCellExperiment that plotPCA could use, at least it is not called PCA which by default is expected. Use runPCA first, please read the manuals of scater and scran first, they have step-by-step code.


The various Gene Ontologies are arranged as Direct Acyclical Graphs. This means that each term can have children and parents. When a gene is annotated with a GO term, it also inherits all of that term's parents. Therefore, when you run a GO enrichment analysis, it is very common that you find the broader terms enriched since those are the ones shared by most genes.

For example, the Biological Process (BP) GO term GO:0003163 sinoatrial node development has the following tree:

￼ GO:0008150 biological_process
     ￼ GO:0048856 anatomical structure development
     ￼ GO:0032501 multicellular organismal process
         ￼ GO:0032502 developmental process
         ￼ GO:0007275 multicellular organism development
             ￼ GO:0048731 system development
             ￼ GO:0009888 tissue development
                 ￼ GO:0048513 animal organ development
                 ￼ GO:0072359 circulatory system development
                 ￼ GO:0060537 muscle tissue development
                     ￼ GO:0048738 cardiac muscle tissue development
                     ￼ GO:0007507 heart development
                     ￼ GO:0014706 striated muscle tissue development
                         ￼ GO:0003161 cardiac conduction system development
                             ￼ GO:0003163 sinoatrial node development
                                 ￼ GO:0060921 sinoatrial node cell differentiation

This means that any gene annotated with GO:0003163 will also be annotated with every one of the terms above it in the tree, all the way to the top of the graph. All genes with at least one BP GO term will also be annotated with GO:0008150, the root of the tree.

As you can see, your "muscle tissue development (GO:0060537)" is also there and is a parent of "sinoatrial node development". However, if you just ran an enrichment analysis (I ave no idea what you did since you haven't actually told us anything), then what will have happened is that there are many genes annotated with various different child terms of GO:0060537, but since all child terms will inherit the parent annotations, what you see is the parent and not the specific term you are looking for.

So yes, there are far more specific terms that what you show, but whether or not they appear in your results will depend on what you are doing.


How to link GDC ids to CCLE cell line names?

The filenames of those bams contain the cell line names. They are constructed:

SAMPLEID.CELL_LINE_NAME.REPLICATE.bam

Thus the first 5 cell lines in your screen capture are SK-N-F1, TE-9, CJM, HCI-H1915 and CAL-33.

All the information that is need is also present in the experimental metadata xml file attached to each sample record. This is accessible by downloading it from the record on the website, or by querying the GDC metadata API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api).


To add a more complete answer: the current coronavirus is closely related to the SARS virus that caused the outbreak in 2004, and on which much research has been done.

Here is a general review of the coronavirus epidemiology, life cycle etc.

I haven't found yet any materials about the RNA structures in the translatable region, however the structures in the 3' and 5' terminal regions are well studied and are known to be crucial for initiating viral replication and for its packaging.

This recent article addresses in details the RNA structures in 5' terminal of human coronaviruses, whereas this older article describes the structures both in 3' and 5' terminals.

Finally, the interesting algorithmic developments for the link that I added to the question are described here.


Declaring an array of strings and integers

When you read in your two CSV files, you can explicitly declare the variable types for each column:

> a <- read.csv(..., colClasses=c("character"), ...)
> b <- read.csv(..., colClasses=c("character", "numeric", "numeric", ...), ...)

You could also investigate as.is or stringsAsFactors, although this can cause problems where there are duplicate strings. Using colClasses is more explicit.

Read the help documentation in ?read.csv for a more complete explanation of these operators.

You can filter with the %in% operator in match:

> strings_of_interest <- as.vector(a$V1)
> b.filtered <- b[b$V1 %in% strings_of_interest,]

If I am reading your question correctly, perhaps you can then take b.filtered and turn that into the matrix you want.

Avoid use of for loops where you can. They can slow R down considerably. A good read on this subject is in the third chapter of Patrick Burns' The R Inferno.


So consider this, you know how many reads map to each of your genes (1). Does that tell you anything about (2) the number of reads that mapped against non-gene regions (3) the number of reads that did not map against the genome at all?

If you only use the gtf files, you know the number of (1) but not of (2) and (3). To me it sounds like you would want to do something as described in this post.

What it boils down to is:

Finding out how many (unique) reads are mapped, counting pairs once

samtools view -F 0x4 foo.sorted.bam | cut -f 1 | sort | uniq | wc -l

Counting how many total reads you have. Since you have paired-end data the number of forward reads suffice:

cat forward_reads.fastq | grep '@' | wc -l

The percentage is then simply mapped_reads/total_reads


Use of heterozygous SNPs in cancer research: why?

I’m no longer working in tumour sequencing so I’m by no means an expert. But in a nutshell, the reason is that, as indicated, homozygous SNPs aren’t informative: if your allelic fraction is 100%, meaning you only observe a single nucleotide at a given position, we don’t know whether we’re dealing with a tumour related mutation. In fact, we probably don’t because the sample likely contains healthy as well as tumour cells.

By contrast, tumour-specific variants always have an AF < 100%, in other words they exhibit some degree of heterozygosity.

In fact, in healthy tissues you’d expect to see one of three cases: your variant allele (= sequence base ≠ reference) has a frequency of either 0% or 100% (homozygous), or 50% (heterozygous), ignoring rare somatic mutations. With tumour samples, by contrast, the fraction will vary depending on the penetrance of the variant in the tumour, and on the fraction of healthy to tumour cells in the sample.


I would recommend looking into Python's plotting utilities, for example here.

I would recommend specifically looking at the examples that plot mathematical functions.

The basic idea is that you create a function that computes your equation based on the different parameters as inputs ($a, a_{max}, l_{max}$) and then you can plug in a range of values for each parameter and plot e.g. with $a_{max}$ on the x axis, then different lines showing a range of values for $l_{max}$, or whatever makes sense for you.

e.g.

import numpy as np
import matplotlib.pyplot as plt

def fn(a, amax, lmax):
  return lmax/2 + (a/amax) * (lmax/2)

plt.figure()
amax = np.arange(0.01, 1, .01)
plt.plot(amax, fn(.005, amax, .1), "ro") ## red points
plt.plot(amax, fn(.005, amax, .2), "go") ## green point
plt.show()

I have no idea if these parameter ranges are meaningful for you, you will have to figure that out for yourself.


full visualisation of draft genomes alignment

Found a solution, using D-Genies, worked great.

Some examples from their website:

Thanks to @user172818.


I am the creator of PixelLib. PixelLib supports instance segmentation of 80 classes of objects using pre-trained coco model. Cell organelles is not part of the classes supported.
The only solution is to train your own custom model. PixelLib supports custom training of a model using your own labelled dataset. Read this article I published on how to train a custom model with PixelLib.


How can I call structural variants (SVs) from pair-end short read resequencing data?

I think the best method or combination of methods will depend on aspects of the data that might vary from one dataset to another. E.g. the type, size, and frequency of structural variants, the number SNVs, the quality of the reference, contaminants or other issues (e.g. read quality, sequencing errors) etc.

For that reason, I'd take two approaches:

Try a lot of methods, and look at their overlap
Validate a subset of calls from different methods by wet lab experiments - in the end this is the only real way of knowing the accuracy for a particular case.


less than 80% of your list has mapped to your chosen identifier type

DAVID last received an update in May 2016 (https://david.ncifcrf.gov/content.jsp?file=update.html), so any additions to the Gene Ontology database since then won't be present. I checked the examples you posted in Biomart and found that only 2 have GO annotations, so you're likely not missing an immense number of genes. Having said that, you'd be best served to use a more frequently updated tool than DAVID. Around 2015 there was an outcry because it hadn't been updated in 6 years and the database was so out of date that people were getting questionable results. Have a look at one of the many Bioconductor packages for convenient free solutions that are better maintained.


Use GSEMatrix=FALSE, the information you want appears to be in the soft file, rather than the matrix. Yes, this is moderately annoying since it means you need two objects to describe this particular dataset.


The only bioinformatics tool other than GATK4 that I am aware of that uses spark is Hail (a Spark based replacement for Plink). Hail is also supported by researchers at the Broad.

Most places I have worked at have not switched over to Spark.  As such, I don't think it is widely used generally. Therefore, I don't think knowing Spark qualifies as a necessary skill to learn for bioinformatics at this time.

However, I think Spark is a superior paradigm for handling the large amount of data that bioinformaticians routinely deal with, and I think we will soon see the field move towards using it more.  Any bioinformatician would do well to acquaint themselves with Spark and to play around with Hail.


factoextra: Error in if (xlab == "x") xlab <- "x value" : argument is of length zero

I actually found an answer just accidentally. It is very unfortunate that factoextra documentation does not explicitly say that data parameter should be of data.frame type. In my case rtsne$Y was a matrix. The error that factoextra is giving is very uninformative.

Here is the code that ultimately worked:

fviz_cluster(db,
----->  data = as.data.frame(rtsne$Y),
        stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
        geom = "point",palette = "jco", ggtheme = theme_classic())


Folded Protein Chunk Dimensional Classification?

CATH
The CATH database classifies protein by fold: https://www.cathdb.info/
So the value from that is probably the most useful for you.
Crystallographic conditions

pH
Salinity
Dissolved Sugars
Reducing Agent Concentration

Your crystallising conditions do not mean much. They are a solution that is close to precipitating your protein but slowly enough for it to crystallise. A protein can crystallise under multiple different conditions and a crystal may be chosen simply because it is a lower index well in the Hampton screening conditions...
These mean very very little. And unless you are planning on doing ML to determine the best screening condition these should be avoided.
Data from the PDB
You have probably noticed that the only stats PDB/PDBe give are crystallographic, like space group, unit cell dimensions, resolution, completeness in shells, Rfree, Rmerge etc. So nothing there is of much use... except for bound ligand. Often this is the cofactor, but often it's a cofactor analog or a product (as opposed to the substrate). But that information can come from Uniprot etc.
pH
A solution has a pH, a protein has a pI. At a pH equal to the pI it will be neutral —and less insoluble. ProtParam can calculate the pI from a linear sequence. However, surface charge totally shows different spots. APBS can map the solution to the Poisson-Boltzmann equation to your protein. There is a PyMOL wrapper, a Pyrosetta wrapper and more. But you still would need to some stats and 3D operations to convert it into a small vector of values —rototranslations and projects are not for the faint of heart or in a rush.
Melting temperature
BRENDA has lots of enzymatic data on enzymes including melting temperature or kinetics at different temperatures (which follow the Eyring equation (TST), until MMRT kicks in) which could be useful.
For a given protein you can calculate the relative folding ∆∆G (eg. with Rosetta Score), which is a proxy to melting temperature in a buffer it is happy in.
Despite paining me to say... some fanciful Poisson-Boltzmann map derived vector and ∆∆G are probably all you could get from analysing structures. Whereas the rest are from CATH classification and other databases that are not structure based.


Most genomics data that is published should be available in GEO (or SRA). Usually in the paper they mention the accession number. However, not every study ends up in there, and sometime studies that are in GEO are not reproducible (poor description of the samples, or design). But if you look at the original papers, most of the times you'll find an accession number.


Salmon output: transcripts quantified with zero reads support

All transcripts should always be returned in the output, whether there's evidence supporting their expression or not. This is for the sake of simple convenience, since typically one would be merging these values across samples into a matrix. Missing values will then tend to become either NA or cause an error, when what would be appropriate would be to have values of 0. This also facilitates calculations such as "what's the median number of reads per transcript?" or "what percent of transcripts have an expected reads count of at least X".


How to concatenate "by chromosome"-VCFs?

I would recommend bcftools concat.  You can't just cat them together because each file has a header section.  The bcftools command will handle all that for you.  Each vcf file must be sorted prior to calling concat

bcftools concat -o total_chroms.vcf chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf


Efficiently concatenating disjoint, ordered bigBed files

No, there's no simple way to do that without someone writing a program tailored toward it that handles the footer section and recreates the zoom levels that follow the compressed blocks of intervals. If you wanted to write such a program, modifying bigWigCat from Kent's tools would be a reasonable tack.


What is the difference between the entity_src_genCategory tag and the entity_src_natCategory tag in pdb-xmls?

src_nat = isolated directly from a natural tissue

src_gen = the source was genetically manipulated

The files that you use (PDBML/XML) are derived from the primary PDB data format which is PDBx/mmCIF.
The names in XML seem to correspond to mmCIF names, so you can read documentation on http://mmcif.wwpdb.org/, and here are direct links:

category entity_src_nat
category entity_src_gen

Can they occur together in one file? I suppose they can for complexes: one "entity" could be natural and another genetically manipulated.


alternatives to MEDIPS to analyse MeDIP datasets

I guess QSEA should be an obvious answer:

  QSEA (quantitative sequencing enrichment analysis) was developed as
  the successor of the MEDIPS package for analyzing data derived from
  methylated DNA immunoprecipitation (MeDIP) experiments followed by
  sequencing (MeDIP-seq).

See http://bioconductor.org/packages/release/bioc/vignettes/qsea/inst/doc/qsea_tutorial.html


Is it possible to correct bond order and conectivity problems in an sdf file according to mol2 files?

RDKit has a one line command to fix the bond order (docs)
 mol = Chem.MolFromMolFile('wonky.mol')
 template = Chem.MolFromMol2File('template.mol2')
 newMol = AllChem.AssignBondOrdersFromTemplate(template, mol)

The problem is that the two need to be identical. I see a carboxylate in your screenshots, so you most likely did the correct thing and docked with the appropriate protonation for pH 7, i.e. as the conjugate base with a negative formal charge. If one does not then RDKit will fail.
Also, RDKit is "not great" with mol2 files. If you have the original SMILES that is better. If the problem is atom names, using a PDB is a better call. If you want to assign atom names, this blog post of mine may help.


Question regarding the function pairwiseAlignment in R

According to the help page, it seems more logic to put the reference sequence in subject, and your query sequence in pattern, but it won't really matter for the outcome (except that the pattern sequence will be on top of the alignment).

For your second question:

  the scores for the optimal pairwise alignments is returned


storing SNPs in the genotype file for each chromosome in a separate file

I believe a pandas approach would be fast enough (probably faster than any other Python approach without pandas).

After converting your table to a pandas data frame, you can do something like df.filter(regex=(".*Chr01")) for the first chromosome then write the resulting data to file with the to_csv() method.

import pandas as pd

data = pd.read_csv('input.txt', delimiter = "\t", index_col = False)

for i in ["Chr01", "Chr02", "Chr03"]:
    regex_pattern = ".*" + i
    print(data.filter(regex=regex_pattern))


What are the p01sum and p0prop parameters in godon?

Disclaimer: I'm the author of godon. The question above comes from the real communication with a software user.

In the godon implementation of the branch-site model I performed the following reparametrization:

$p_{01sum} = p_0 + p_1$

$p_{0prop} = \frac{p_0}{p_0 + p_1}$.

This has the following advantages:

Both parameters have a well defined range $(0, 1)$.
There is less dependency between the two. E.g., $p_2=1-p_0-p_1$ (proportion of sites under positivie selection) does not depend on $p_{0prop}$.
Both properties are very helpful for both likelihood maximization and MCMC.

In case you are interested in $p_0$ and $p_1$ there is a straightforward way to go back to the original parameters:

$p_0=p_{0prop}*p_{01sum}$

$p_1=p_{01sum}-p_0$.

P.S. I also added this information to the tutorial.


Can I get the graph generated by cellranger

Cell Ranger

You can't download the tSNE coordinates for cells directly from the Analysis tab of the fancy, polished .html document that Cell Ranger produces. If you have access to the machine on which the pipeline was run, you can grab the intermediate files. Everything is well documented by 10X. They explain defaults etc. on that page. The tSNE coordinates are available at:

analysis/tsne/2_components/projection.csv

and the contents look like:

$ head -5 analysis/tsne/2_components/projection.csv
Barcode,TSNE-1,TSNE-2
AAACATACAACGAA-1,-13.5494,1.4674
AAACATACTACGCA-1,-2.7325,-10.6347
AAACCGTGTCTCGC-1,12.9590,-1.6369
AAACGCACAACCAC-1,-9.3585,-6.7300

Seurat

And the second point: Cell Ranger is replicating the Seurat workflow. Seurat uses the Louvain method for community detection by default. Check out the intro tutorial and see if you can follow along. This tutorial works with 10X data, in fact the first function call in the tutorial reads in the 10X matrices. If you are already plotting in R you'll speed through this tutorial. Seurat is great to work with, you'll have tSNEs and clusters in no time. And there's so much more you can do from there!


Recently I've developed a python lib for converting InChI to InChIKey: https://github.com/liwt31/chembl_ikey. Maybe it will help.


AttributeError: module 'scater' has no attribute 'normalize'

I was able to get to the normalized values using logcounts function in the following way:

mtx_norm_r = sceLib.logcounts(sce)
mtx_norm = np.array(mtx_norm_r)

Since I used log_exprs_offset=1 the field that I needed to use is logcounts, otherwise it should be normcounts which I found here.


How is a principal component analysis for population structure made from a bam or other alignment file?

In general, the process works as such:

Variants are found at a number of loci in each sample (missing loci are often imputed).
Each sample is scored at each locus, often according to the number of copies of the reference allele that it has (so, 0-2).
The resulting matrix is used for PCA.


how to repeat records in fastq n times efficiently?

Not sure why you would want to do this, but if the aim is to repeat all records in the fastq/fasta files n times, the most efficient way I can think of is just to concatenate the orginal file n times, similar to this thread.
i.e adapting the accepted answer would look something like this (I have not tested this):
for i in {1..n};do cat file.fq >> file_ntimes.fq ; done

Update based on comments:
To replicate each fastq entry back to back, maybe this hacky unix command can be adjusted (say for n=4):
n=10
paste - - - - <file.fq | awk '{for(i=0;i<$n;i++)print}' | tr '\t' '\n' > file_10times.fq

The idea is to paste blocks of 4 lines (i.e. 1 fastq entry) into one line and replicate this line n times. tr is reversing the paste action (since the separation by paste is a tab character). This won't work if fastq file for some reason has tabs in it which I don't expect but didn't test either.
Note, for fasta file this would be paste - - rather than paste - - - -
I haven't tested this, but if it works for you I could polish it up if needed.
Update2
Based on Xioafeis comment this is solution seems to work as well:
bioawk -cfastx '{for(i=1; i<=10; i++) {print "@"$name"\n"$seq"\n+\n"$qual}}' file.fq > file_10times.fq


What is the best way to account for GC-content shift while constructing nucleotide-based phylogenetic tree?

There are models that take into account compositional heterogeneity both under the maximum likelihood and Bayesian frameworks. Although the substitution process is not time-reversible, the computations are simplified by assuming that the instantaneous rate matrix can be decomposed into an "equilibrium frequency vector" (non-homogeneous) and a symmetric, constant exchange rate matrix.

I guess all your suggestions are also valid, and I remember recoding being used successfully to reduce the GC-content bias (examples in the references above and here).


How does Li and Durbin's BWA paper compare alignment programs on real data?

What they are pairing is the RNAseq data with the reference genome: as the legend of the table 2 says:

  The 12.2 million read pairs were mapped to the human genome

These softwares align the same read to the same genome. They aren't comparing if it is mapped to the same location, only how many mismatches has each program.


What are some good practices to follow during EPIC DNA methylation data analysis?

EPIC data can be processed in the same manner as the previous iteration of methylation array data from Illumina (450k). This means that starting with .idat files, normalization should be performed (for example, via the minfi package). A recent paper from the creators of minfi is particularly helpful because it makes clear that normalized EPIC data from their package can be immediately compared against, for example, level 3 TCGA data.

After that, I suggest using the manifest to attach genomic coordinates to your probes and segregating them into functional regions. By testing differential methylation in only regulatory regions, for example, you can increase the statistical power by reducing the overall number of tests to the ones you expect to yield major differences.

There are existing packages out there for differential methylation analysis, but without knowing your replicate structure or aims, it is difficult to point you in the right direction.


Assuming that gap filling is a domain-independent and works the same for viruses as for bacteria or eukarytes.

There are two real options I am aware of: Sealer a successor of GapFiller from ABySS package and GapCloser from SOAPdenovo2 package. According to Sealer paper it's better, but I am not sure if any indepented benchmarking have confirmed their conclusion. On my insect data GapCloser worked better than GapFiller (back then I have not tried to compare it to Sealer).


Biopython BioSQL error : AttributeError: 'DBSeq' object has no attribute '_data'

It's a bug in Biopython 1.78.
What those three sequences have in common is that they're less than 61 characters long.  When you print the record it implicitly calls __str__ for the record, which calls __repr__ for the Seq.  DBSeq gets its method from Seq, which looks like:
    def __repr__(self):
        """Return (truncated) representation of the sequence for debugging."""
        if len(self) > 60:
            # Shows the last three letters as it is often useful to see if
            # there is a stop codon at the end of a sequence.
            # Note total length is 54+3+3=60
            return f"{self.__class__.__name__}('{str(self[:54])}...{str(self[-3:])}')"
        else:
            return f"{self.__class__.__name__}({self._data!r})"

...and then it only tries to access the nonexistent _data for the short sequences.
So there's a mismatch between DBSeq's lack of a _data attribute and its use of Seq.__repr__.  I don't see a DBSeq class at all in the latest Biopython, so maybe this is fixed already.  But there's no tag for 1.79 yet that I can see and I can't get the very latest code to import successfully so I'm not sure.  In the meantime if you avoid letting repr(record) happen you should be able to work around it.


Thanks to the helpful comments I figured it out:

samtools view -b -h -L bedfile.bed originalbam.bam > newbam.bam

-b = Output as bed file
-h = include header
-L = only output alignments overlapping the input bed file


How does the BWA-MEM algorithm assign its mapping qualities?

Yes, there bwa-mem was published as a preprint

  BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference
  position x with the best extension score achieved at query position y.

  ...

  Secondly, while extending a seed, BWA-MEM tries to keep track of the
  best extension score reaching the end of the query sequence

And there is a description of the scoring algorithm directly in the source code of bwa-mem (lines 22 - 44), but maybe the only solution is really to go though the source code.


Identifying 15 clusters from 500 cells is challenging. It will depend on the relative abundance of each cell type within the tissue. You can use the cell type frequencies found in the shortSeq data set to see what to expect from a 500 cell dataset using the howmanycells tool

About the pc plots:
The first is an elbowplot which explains how many variation is explained by each PC. The knee basically defines when components become less informative.

EDIT: I noticed that the y-axis of the elbowplot starts at 4 suggesting that the dimensionality of your data is most likely ~10 dimensions. But I expect many of the clusters that don't show up on your Umap plot to be present at very low frequencies.

The second plot is a jackstraw plot and provides a visualization of a statistical test that calculates the significance of your principal components. This one shows that all your components are significant. You can combine this with the elbowplot to estimate the dimensionality (10-15 dimensions) of your data and decide which components to use for downstream analysis. In general the more they dropoff to the dotted line the less significant they become.

While these plots are useful as a first impression, always rely on biological information and common sense to decide what information (here: components) to include in downstream analysis such as Umap, clustering, etc. For example: usually you don't want to include components mainly defined by ribosomal, mitochondrial or hsp genes.


Software for microbial profiling from 16S rRNA gene sequence

I tried dada2 and is not bad (if you know R).

QIIME2 is also an option.
Many other are available, the choice might also depend on your sample and your exact question.

For functional profiling you may try picrust2


How to map multiple CDS sequences to exome data using bwa mem?

If I have one file containing all the CDS from my reference species, can that be used in bwa-mem?

That's fine, normal usage of bwa-mem entails giving it a fasta file with all chromosomes/contigs/transcript/etc. in an organism, so your usage is no different. How your custom "make an MSA" script handles that is unknown of course.


You're supplying a GTF feature file, but telling it that the features are in GFF3 format. I think you need -fformat gff2. See https://www.ebi.ac.uk/Tools/sfc/emboss_seqret/help/


Use grep to get SNPs that have ";" in the name:

>$ grep ";" myMap.txt | cut -f1 -d" " > myBadSNPs.txt

Then use plink with --exclude :

plink --file myFile
--exclude myBadSNPs.txt
--out myFileFiltered


Try with:

makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90

The -out is the name of the DB. So if you use uniref90.fasta, it creates uniref90.fasta.pinwhile the script seems to be expectinguniref90.pin`.


Do you know a program/script to count how many sequences have mismatches on each position of my primers?

I'm amazed at your tenacity in getting Excel to do what you want; that looks like many hours of work to count mismatches.
The hammer I use at the moment for this sort of stuff is LAST, which will do local alignments of sequences to references. If the number of mismatches are not too large (as would be expected for primer sequences), a local alignment will extend to the end of one of the sequences, and I expect that should work fine for your purposes.
I'm going to demonstrate two options, one that can be done just using Excel and web services, and one that needs command-line access. Note that these are only two ways to fetch mismatch counts; you might want to do more than that (it's not clear from your question), and there are many other ways. Also, this may not work for your particular circumstances, but I think it should be quick enough that it's worth a try:
Option 1 - Web-based LAST
In this, I use the Web version of LAST to generate alignments for the primers [as "reference"] mapped to the sequences [as "query"].

To help create the mismatch (or in this case more correctly, "edit distance") statistic, in the Advanced Parameters section, set the a: gap existence cost to 1, and the b: gap extension cost to 1.
Click on Align!
Once the Alignment Results appear, click on the Alignments in tabular format link down the bottom, and save the file to your computer as last.tab.csv (the .csv at the end is the most important bit, to make sure that Excel recognises it as something it can open).
Open that file in Excel, set the delimiter to <tab>
The first column should be the alignment score, which in this case is the match score. This can be converted into a mismatch score by subtracting it from the primer sequence length (column 6).

Option 2 - Local LAST Installation
Installation
Install LAST and samtools.
Index Generation
The first step is to generate a reference index. LAST doesn't care too much about what's query and what's reference, but because I want all the genome sequences to appear on each line, the genome sequences will be the query, and the primer sequences will be the target. I use Martin Frith's lastdb mantra to do this:
$ lastdb -uNEAR -R01 primers.fa primers.fa

Mapping
The next step is to do the actual alignment with LAST. I set the match and mismatch to be both 1 (i.e. the default for the Web version), set the gap existence and extension to be 1, change the output format to tab, and save the result to last.tab.csv:
$ lastal -r 1 -q 1 -a 1 -b 1 primers.fa sequences.fa -f tab > last.tab.csv

Aggregating data
Finally, it's a matter of aggregating the count information. I use grep -v to exclude lines that begin with a #, use awk to subtract field 1 from field 6, then sort the results and use uniq -c to count the mismatches:
$ grep -v '^#' last.tab.csv | awk '{print ($6 - $1)}' | sort | uniq -c
2118 0
   5 1
   5 2


How can I specify to DEseq2 to only perform comparisons on which I am interested with?

You can specify the exact comparisons you want in the results() function. So:

dds = DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds = DESeq(dds)
res = results(dds, contrast("condition", "treatment1", "control1"))

Note that condition should not be, "stim001-control1, stim001-control2, etc.", but instead, "stim001control, stim001control, stim001control, stim001treatment, etc.". Don't put minus signs in your levels and don't number your samples in them (they're not in different groups.

The last command would be repeated for each of the comparisons. Note that this will still be quite slow due to the number of samples and the number of groups. With so many samples, you might just make a separate sampleTable for each of the comparisons you want to make. Each would then only contain the samples relevant for that comparison, so 60 total in each.

Alternatively, if you want to fit all of the groups at once then it's likely that limma/voom will prove to have better performance. It uses different math that happens to be quicker with large models.

BTW, if the slowness is coming from DESeqDataSetFromHTSeqCount() then you can mere everything on the command line into a single matrix. A couple lines of python or even the join command can do that quickly enough.


Read sorted and indexed BAM files faster in C++?

In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest.

In the first method I use an std::istreambuf_iterator, in the second I read the file line by line into a std::vector and in the third one I use C flavoured operations.

/*
compile with g++ -std=c++11 -O2 reading_files.cpp -o reading_files
test the program with the following
dd if=/dev/random of=./dummy5GiB bs=1024 count=$[1024*1024*5]
 for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
*/

#include <fstream>
#include <iostream>
#include <vector>
#include <cstring> //for strcmp
#include <string>
#include <algorithm> //for min

using namespace std;

int main(int argc, char * argv[])
{
    if(argc != 3) {
        cerr << "usage: reading_files <path> <number 0 or 1 or 2>\n";
        return 1;
    }

    ifstream file(argv[1], ios::in|ios::binary);
    if(!file.is_open()) {
        cerr << argv[1] << " file not opened!\n";
        return 1;
    }
    cout << "---------------------------\n";
    if(!strcmp(argv[2],"0")) {
        cout << "method 0\n";
        vector<char> content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
        cout << content.size() << " bytes read\n";
    } else if(!strcmp(argv[2],"1")) {
        cout << "method 1\n";
        string line;
        vector<string> content;
        while(getline(file,line)) {
            content.push_back(line);
        }
        cout << content.size() << " lines read\n";
    } else if(!strcmp(argv[2],"2")) {
        cout << "method 2\n";
        file.seekg(0, std::ios::end);
        unsigned long long file_size = file.tellg();
        cout << "file size is " << file_size << "\n";
        //maximum 5 GiB of memory, so to fit my RAM
        unsigned long long max_block_size = 1024ULL * 1024ULL * 1024ULL * 5;
        file.clear();
        file.seekg(0, std::ios::beg);

        unsigned long long block_size = min(max_block_size,file_size);
        char * buffer = (char*)malloc(block_size*sizeof(char));
        unsigned long long processed = 0;

        while(!file.eof()) {
            file.read(buffer, block_size);
            cout << "read a block of " << min(file_size - processed, block_size) << " bytes\n";
            //do work here with the block in memory (if the block size is > 0)
            processed += block_size;
        }

        free(buffer);
    }
    file.close();
    return 0;
}

In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.

Here is the benchmark of the program on a 5 GiB randomly generated file.

bash-3.2$ for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
---------------------------
method 0
5368709120 bytes read

real    0m35.194s
user    0m19.665s
sys 0m12.999s
---------------------------
method 1
20973135 lines read

real    0m40.478s
user    0m34.046s
sys 0m5.870s
---------------------------
method 2
file size is 5368709120
read a block of 5368709120 bytes
read a block of 0 bytes

real    0m4.413s
user    0m1.757s
sys 0m2.617s

Here is the benchmark on a real 32GiB .bam file using only the third method, here the splitting of the file read plays a fundamental role. I can't tell you how big because my hard disk memory gets saturated before the end of the execution and so I have to kill the program.

time ./reading_files ./HG00252.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam 2
---------------------------
method 2
file size is 34316054058
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 2103799338 bytes

real    0m50.656s
user    0m9.135s
sys 0m35.728s

The key point is to minimise the file operations reading large block of data at once (as long as they fit the memory) and operate on them.

Using mmap could lead to an improvement, but I don't think it would be as large as using the third method instead of the first two.


How can I restrict the search space when querying GEO?

You can use the following query to first search GEO to get a list of GEO Series of interest, then find all linked SRA runs and perform a second filter step to keep only the SRA runs that satisfy another set of criteria as shown below.
esearch \
  -db gds \
  -query "daf-2[all fields] \
    AND expression profiling by high throughput sequencing [DataSet Type] \
    AND Caenorhabditis elegans [organism]" \
  | elink \
  -db gds -target sra \
  | esearch \
  -query "(#3) \
    AND Caenorhabditis elegans [organism]" \
  | efetch -format runinfo

Here, I am combining two independent queries using the (#3) term in the second query. More information about this is in the "Combining Independent Queries" subsection of "Searching and Filtering" section of the Entrez Direct documentation.


Theoretical limit of human genome compression

TLDR; For the most modern CHM13 assembly (v1.1; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best compression without any use of a reference would be between 624 MB and 641 MB.
Upper estimates of genome compressibility can be made using standard compressors (see, e.g., this comparison). In my experience, xz can achieve great compression ratios on genomic sequences, especially when used with the maximum compression level (-9) and a single thread (-T1). Importantly, the original FASTA file needs to be cleaned – linebreaks hamper compression as they act as noise; the reformatting can be done by seqtk seq -U (which also converts all nucleotides to the upper case). It the case of the CHM13 assembly, xz compresses the human genome to 641,369,600 bytes.
Lower estimates are trickier. The most straightforward approach is through the number of distinct k-mers. Since k-mer sets contain less information than the original assemblies, we can use them for determining the lower bounds of genome compressibility. The lower bound on the compressibility of k-mer sets is 2 bits per distinct k-mer (under certain assumptions that I believe are satisfied in our case). The CHM13 human genome assembly contains 2,497,988,809 distinct canonical 31-mers, which corresponds to 624,497,203 bytes.
Overall, xz does a great job at compressing single human genomes, and the resulting size is very close to the optimum. The full analysis from this post can be found on http://github.com/karel-brinda/human-genome-compression. The resulting files are here: https://zenodo.org/record/4776580.


pysam or piping samtools view to a python script

It's a bit hard to say with certainty, though I would suspect that offloading the BAM decompression by using a pipe will be very slightly faster. Note that decompressing and parsing the BAM file will not be the bottleneck in your processing, rather the python script itself will be.


Genotyping technologies do not maintain phase?

Phasing means that you are able to distinguish the both copies of your target sequence. (Usually you have one from your mother and one from your father).

If you can do this, depends mainly on the sequencing technology you are using. For example in case of sanger sequencing, you cannot do this, because you (pre)amplify all your copies at once and detect the sequence at the same time at the same location.

In case of ngs you have read information at the end. One read represents the sequence of one original dna molecule. If you detect multiple variants in one read, you can be sure that they are from the same chromosome and thus they are in phase. Now you can go on and try to find other reads that have unambiguous overlaps with these variants on your read. If you find some, you can also say that these variants have their origin on the same chromosome.


An ent file is a pdb file. Just change the extension. PDBe for example provides their PDB files as such.


What is a simple command line tool for doing Needleman-Wunsch pair-wise alignment on the command line

You are looking for the needle program from the EMBOSS suite. Available in bioconda.

http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html

To read sequences from the commandline, you need so specify the format as asis. To get the output on the screen, you'll need -stdout and to use the default alignment parameters (gap penalty 10, extend_penalty 0.5) you'll need -auto. Thus your query above would be:

$ needle -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG -auto -stdout
########################################
# Program: needle
# Rundate: Tue 27 Nov 2018 10:48:50
# Commandline: needle
#    -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
#    -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG
#    -auto
#    -stdout
# Align_format: srspair
# Report_file: stdout
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: asis
# 2: asis
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 56
# Identity:      37/56 (66.1%)
# Similarity:    37/56 (66.1%)
# Gaps:          18/56 (32.1%)
# Score: 181.0
#
#
#=======================================

asis               1 GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---     47
                              |||||||||||||||||||||||||||||||||||||.
asis               1 ---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAG     41

asis              47 ------     47

asis              42 AGGAGG     47

#---------------------------------------
#---------------------------------------


I am getting this error when running a test data set for Lyve-SET on RHEL 7. Is this an MPileup error or something else?

The samtools v0.1.19 manual confirms that there is no --region option.  In that version there is only the -r argument to pass in the targeted region.  If you update to a newer version of samtools this should execute correctly


Evidence of emergence of "genuine" novel protein interactions by aa mutation

HA Sialic acid receptor binding in influenza virus determining host specificity, viz. Bird flu vs human flu.


The expression values for each gene are scaled / standardized by subtracting the genes mean expression and dividing by its standard deviation. A value of -1 would imply it's one standard deviation below the mean expression for that gene.


Get gene sequence based on the annotation

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
```


Calculate the percentage of each unique phylogenetic tree in a BEAST output

I finally managed to do it in R. Here is my code:

install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)


You don't have a RAM problem because you only have about:

  10-20 reads that I want to access from each file.

Please go for (2). If you are concerned about RAM usage, you may selectively write some of your data back to disk.

Tens of thousands of 10-20 reads is not a problem for modern computing.


Getting hierarchy of cell populations with Drop-seq data

This article uses the freely available R package dropbead for filtering and then Seurat to perform a principal component analysis that groups together affine transcriptomes. It could be what you are looking for.


Identifying mutually **exclusive** gene sets

"Mutually exclusive" is not a precise statistical term because it could be seen as independence.

Anyway if truly mutually exclusive AND their behaviour follows periodicity, i.e. sine waves, which it might do in a cell cycle, i.e. the cells divide at periodic intervals, you can do this via asynchronicity.

If the gene expression can be described under periodicity analysis you simply sum the amplitude between two genes. When the amplitude is zero ... they are expression is antogonistic - when one is expressed the other is silent.

Amplitude simply means the height of the period.

Periodicity analysis is quite simple if you have uniform time measurements.

If you have 'missing data', i.e. you time intervals are not uniform, for example you skipped a reading then this is called missing data in periodicity analysis. Missing data periodicity analysis is complicated, but its doable.

So in this context its quite simple, you would identify the periodicity of all genes and construct a 2x2 matrix subtracting their respective amplitudes. It would be kinda juicy to do.

The above is cool if the genes are periodic, however covariance is a better solution in general terms.

What you want is negative co-variance. If two genes had a covariance of -1 ... that will mutually exclusive. As above 2x2 matrix and the co-variance across the grid. You do need to be careful that the program you are using doesn't remove negative values of the covariance.

Covariance = 1, meaning two genes are co-expressed
Covariance = 0, random distribution
Covariance = -1 mutually exclusive expression


What should be the numbers in the fastSTRUCTURE file?

Numbers, usually in the form of 0/1/2, is a common way many formats represent bases.
For example, if a SNPs has bases which are A/T, then these are encoded as 0/1, with the reference allele/base taking the 0 and the alternate allele/base taking 1.
However, if you read the documentation, it says:

Specifically, rows in the data file correspond to samples, with two
rows per sample (note that only diploids are handled by this
software), and columns correspond to SNPs. The first 6 columns of the
file will be ignored; these typically would include IDs, metadata,
etc. This software only handles bi-allelic loci. The two alleles at
each locus can be encoded as desired; however, missing data should be
encoded as `-9'.

So you should be able to use your base letters here and it will work OK as long as the columns and rows are formatted correctly.


GRCh38 vcf file with common cancer mutations

Your question doesn't give enough information for a specific answer but this should do for a start.

Get the VCF file describing all variants in Clinvar from NCBI:

wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz

Extract any variants whose VCF line contains the word "cancer":

zgrep -iE '^#|CLNDBN=[^;]*cancer' clinvar.vcf.gz > cancer.vcf

Now, you should probably filter that further based on variant frequency, type of cancer, etc etc but that should be a good start.


Why does this naive Python paired Fastq reader work?

Your code only runs prior to Python 3.7, now it raises RuntimeError: generator raised StopIteration. This is essentially a bug fix described in PEP 479 -- Change StopIteration handling inside generators. You can enable the new behaviour from Python 3.5 via from __future__ import generator_stop.

Note, it's not related to file-handles per se, you can reproduce it with a smaller example with a simple iterator:

it = iter([1,2])

def f(it):
    while True:
        yield next(it)

print(list(f(it)))

You should not write code relying on the bug


Which identifier do I have to use when I want to add a ##Fasta section to a gff?

Unless I'm mistaken, the only variant of GFF that explicitly permits Fasta entries is GFF3, and the data you've shown above is definitely not valid GFF3. The mRNA entry seems to correctly use the = symbol to separate keys from values in each attribute, but the remaining columns all use spaces and quotation marks (GTF style) instead. This is likely to cause all kinds of problems with software packages that take GFF data as input.

But to more directly address your question, there is no "right" way to connect feature entries to sequence entries in a GFF3 file. There are some conventions, but these are all very specific to a small set of tools and community of practice, and it is rare that these conventions are objectively enforced. Some options include:

Placing genomic nucleotide sequences in the ##FASTA section, and then using feature coordinates to pull out mature mRNA sequences and translate them in silico to derive the corresponding peptide sequence. In this case, the Fasta entries would be genomic sequence identifiers that correspond to the first field/column in each feature entry.
Placing mature mRNA sequences in the ##FASTA section, and then directly translating this sequence to derive the corresponding peptide sequence. In this case, you'll want the sequence IDs in the ##FASTA section to match the ID attribute of the corresponding mRNA feature.
Placing the peptide sequences directly in the ##FASTA section. In this case, you would either want to have protein features declared in the GFF3 file (uncommon) or have some sort of protein_id attribute associated with each mRNA feature, in which case you'd want the sequence IDs in the ##FASTA section to match this protein_id attribute.

There are probably other ways of doing it, too. Which way is "right"? The GFF3 specification has nothing to say about it. It depends on the tools you are using to do the subsequent processing and analysis. If you're lucky, those tools might even document what they expect in terms of input.

Let me just say, as a post-script, that I think placing sequences and feature annotations in the same file is a horrible idea. I've never liked this about GFF3.


How to colour multiple residues in Pymol?

Note: as I don't have any working version of pymol at the moment,
I am not able to test the solution properly.
However, an easy way is to first create a selection, and then to color it:

select toBecolored, resi 3-8
color yellow, toBecolored

Regarding the L/H chain part, if they are effectively annotated as L and H, something like the following should work:

select toBeColored, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)
color yellow, toBeColored

Combining both commands is probably possible as:

color yellow, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)

Pymol allows synonyms (eg. resi and i.  or chain and c.).

The following pages may be also helpful:

Pymol wiki: help on Selection
SourceForge

And the two following links, if you want to get a bit deeper in proper programming with python for pymol (as they can provide you with inspiration):

Pymol wiki: plugins
Pymol scripts


Two variants associated with same chromosomal position in genotype data for one individual

It's pretty common for SNPs and indels that map to the same position in the genome to be assigned different rsIDs, particularly when the indel spans multiple base pairs in the reference. You can get more info about these particular variants here and here. Or you can see how both of these variants relate to the reference sequence here.

Without info about the reference allele it's pretty confusing (particularly because they appear to be reporting the opposite strand from the database (G vs. C/T), but I would interpret this to mean that both alleles match the reference for rs267607236 (they have the insertion), but do not match the reference for rs6010164 (G-C pair rather than T-A). This is also the most common genotype at this position, being found in almost 80% of people in most populations.


Prediction of prokaryotic origins of replication (ORI)

I have not seen the papers but I know the technique SVM. This is a machine learning approach its pretty good, but has been around this the 1980s.

Its quite simple known ORIs are given with the associated genetic trait, e.g. GC skew, SVM will then regress the GC skew of a known ORI against the GC content of the genome, there will be a training, parameter optimization and test split in the data, typically 80:20:20. De novo data will then be feed into the trained algorithm, in the same way the test data was, and a prediction will be generated. There are weaknesses, the authors have used only one ML technique, typically numerous will be used e.g. SVM, random forest, lasso regression etc... In addition, the prediction accuracy for GC is low (75%), in addition 500 genomes isn't great for ML, you need 1000s.


Since they share the ensembl_ID column, you can merge them, then assign the symbol column to the rownames, then delete the symbol column.

something like:

merged <- merge(mat1, mat2, by = 0)
rownames(merge) <- merged[,'symbol']
merged[,symbol] <- NULL

https://stat.ethz.ch/R-manual/R-devel/library/base/html/merge.html


Get a list of PDB ids from a list of sequence protein number?

A typical Uniprot URL would look like the example here:

https://www.uniprot.org/uniprot/?query=ID:P26378&format=tab&columns=id%2Cdatabase(PDB)

The response will be:

  Entry Cross-reference (PDB)

  P26378    1FXL;1G2E;

We want our query to include the protein ID, the output format should be tab separated and the only contain the columns id and database (PDB).

To get the PDB IDs you would need to split the response, take the last line, split this line after each tab and split it again after each semi-colon.

Using Python's Requests library the whole code would be.

import requests
import time

uniprot_ids = ['P26378', 'O35433', 'Q02910']
url = 'https://www.uniprot.org/uniprot/'

protein_to_pdb = {}
for protein in uniprot_ids:
    params = {
        'format': 'tab',
        'query': 'ID:{}'.format(protein),
        'columns': 'id,database(PDB)'
    }
    contact = ""  # Please set your email address here.
    headers = {'User-Agent': 'Python {}'.format(contact)}
    r = requests.get(url, params=params, headers=headers)

    protein_to_pdb[protein] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')
    protein_to_pdb[protein].pop(-1)
    time.sleep(1)  # be respectful and don't overwhelm the server with requests

print(protein_to_pdb)

To accelerate things, you could store your protein_to_pdb dict in a file and only retrieve the new ones
Checking if the returned ID matches your input ID is always a good idea


Arlequin files not able to converge beyond 2000 steps for some .arp files

What you appear to have is a Bayesian MCMC for 2000 "steps" which isn't a lot and normally the term 'steps' is used for parsimony, so I'm a bit confused.

However, if a Bayesian analysis does not reach convergence, i.e. duplicate runs do not achieve the same results then the data space has not been solved so you get no result.

9 haplotypes does not sound much to me, it is possible there are too many sequences, it is also possible there is simply not enough variation (9 haplotypes ... how many alleles?).

Suggestions

Arlequin should provide a default convergence value which you can lower to force a result
Decrease the number of sequences but maintain the number of haplotypes
Try their example data set

Do you get any associated file such as posterior probability, any output file with a minus sign in front of it?


Comparing two multi-fasta files of the same set of proteins with parser - to find and count mutations after treatment

from Bio import SeqIO
for normal, treated in zip(SeqIO.parse("/data/statistic/normal_samples", "fasta"),
                           SeqIO.parse("/data/statistic/with_treatment", "fasta")):
    ... do stuff...

That's generally how you zip iterators together in python.


Genome databases of microbial eukaryotes

You can find many protist genomes on ensembl

http://protists.ensembl.org/index.html


Quantifying reads mapping to multiple loci

You almost had the correct python code already, you just need to filter out secondary alignments:

def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        if query.is_secondary:
            continue
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts


How to version the code and the data during the analysis?

There is a couple of points to consider here, which I outline below. The goal here should be to find a workflow that is minimally intrusive on top of already using git.

As of yet, there is no ideal workflow that covers all use cases, but what I outline below is the closest I could come to it.

Reproducibility is not just keeping all your data

You have got your raw data that you start your project with.

All other data in your project directory should never just "be there", but have some record of where it comes from. Data processing scripts are great for this, because they already document how you went from your raw to your analytical data, and then the files needed for your analyses.

And those scripts can be versioned, with an appropriate single entry point of processing (e.g. a Makefile that describes how to run your scripts).

This way, the state of all your project files is defined by the raw data, and the version of your processing scripts (and versions of external software, but that's a whole different kind of problem).

What data/code should and should not be versioned

Just as you would not version generated code files, you should not want to version 10k intermediary data files that you produced when performing your analyses. The data that should be versioned is your raw data (at the start of your pipeline), not automatically generated files.

You might want to take snapshots of your project directory, but not keep every version of every file ever produced. This already cuts down your problem by a fair margin.

Approach 1: Actual versioning of data

For your raw or analytical data, Git LFS (and alternatively Git Annex, that you already mention) is designed to solve exactly this problem: add tracking information of files in your Git tree, but do not store the content of those files in the repository (because otherwise it would add the size of a non-diffable file with every change you make).

For your intermediate files, you do the same as you would do with intermediate code files: add them to your .gitignore and do not version them.

This begs a couple of considerations:

Git LFS is a paid service from Github (the free tier is limited to 1 GB of storage/bandwidth per month, which is very little), and it is more expensive than other comparable cloud storage solutions. You could consider paying for the storage at Github or running your own LFS server (there is a reference implementation, but I assume this would still be a substantial effort)
Git Annex is free, but it replaces files by links and hence changes time stamps, which is a problem for e.g. GNU Make based workflows (major drawback for me). Also, fetching of files needs to be done manually or via a commit hook

Approach 2: Versioning code only, syncing data

If your analytical data stays the same for most of your analyses, so the actual need to version it (as opposed to back up and document data provenance, which is essential) may be limited.

The key to get this this working is to put all data files in your .gitignore and ignore all your code files in rsync, with a script in your project root (extensions and directories are an example only):

#!/bin/bash
cd $(dirname $0)
rsync -auvr \
    --exclude "*.r" \
    --include "*.RData" \
    --exclude "dir with huge files that you don't need locally" \
    yourhost:/your/project/path/* .

The advantage here is that you don't need to remember the rsync command you are running. The script itself goes into version control.

This is especially useful if you do your heavy processing on a computing cluster but want to make plots from your result files on your local machine. I argue that you generally don't need bidirectional sync.


If you really just want the start and end positions of all genes in a list of genes from Homo_sapiens.GRCh38.93.gtf.gz, you could do this in awk. For example, if your list of genes is:

$ cat genes
ENSG00000131482
ENSG00000223972
ENSG00000228794
ENSG00000187961

You can do:

$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$4,$5
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 11869 14409
ENSG00000228794 825138 859446
ENSG00000187961 960587 965715
ENSG00000131482 42900797 42913369

Explanation:

if(NR==FNR){ wantedGenes[$1]++ } : NR is the input line number and FNR is the current line number of this input file. The two are only equal while reading the first file. So, if this is the first file, save the first field of each line (the gene name) in the array wantedGenes.
if($3=="gene"){ : if this line is a 'gene' entry.
gsub(/";*/,"",$10); : the gene names in the GTF file are quoted and have an ; after the quote. For example: "ENSG00000187961";. Since we don't have this in the list, we need to remove it so we can match the gene name. This gsub command will remove all " followed by 0 or more ; (the ;* means 0 or more occurrences of ;) from the 10th field of each line. So, it will replace "ENSG00000187961"; with ENSG00000187961.
if($10 in wantedGenes){print $10,$4,$5} : if the 10th field (the modified gene name from above) is in the wantedGenes array, then print the 10th field (the gene name), the 4th (start) and the 5th field (end).

You can easily modify this to print the length instead of the start and end positions:

$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$5-$4
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 2540
ENSG00000228794 34308
ENSG00000187961 5128
ENSG00000131482 12572

As requested in chat, here's a python version of the same thing:

#!/usr/bin/env python3
import re
import sys

wantedGenes = []
with open(sys.argv[1], mode='r') as f:
    wantedGenes = f.readlines()

wantedGenes = [ x.rstrip() for x in wantedGenes]

print('# Gene\tstart\tend\tlength')
with open(sys.argv[2], mode='r') as f:
    for line in f:
        if not line.startswith('#'):
            fields = line.split('\t')
            annotations = fields[8].split(' ')
            geneName = re.sub('[";]', '', annotations[1])
            if geneName in wantedGenes and fields[2] == 'gene':
                print('%s\t%s\t%s\t%d' %
                      (geneName, fields[3],fields[4],
                       int( int(fields[4]) - int(fields[3]) )
                      ))

Save the file as foo.py and run with the file names as arguments:

python foo.py geneList.txt Homo_sapiens.GRCh38.93.gtf


Confirm success or failure of RNA-Seq normalization

You should use box plots and PCA plot. Let's take a look at the RUV paper:

  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4404308/

Before normalization and after UQ normalization:

  Libraries do not cluster as expected according to treatment. ... for UQ-normalized counts. UQ normalization does not lead to better clustering of the samples...

Before normalization, the medians in the box-plot obviously look very different among replicates.

After UQ normalization, the medians look closer but Trt.11 look like an outlier. Furthermore, the treatments aren't clustered on the PCA plot. Since they are replicates, you'd like them be close on the plot.

After RUV normalization

  ... RUVg shrinks the expression measures for Library 11 towards the median across libraries, suggesting robustness against outliers. ... Libraries cluster as expected by treatment. ...

The RUV has made the distribution more robust and the samples closer on the PCA plot. However, it's still not perfect as one of the treatments is not close to the other two on the first PC.

The vignettes for Bioconductor RUVSeq describes the two functions: plotRLE and plotPCA.


How to create Phylogenetic Trees from fasta files in Python or R?

I would not look for a package for this, but instead build a small pipeline calling external tools with something like the following workflow:

Cluster the ~100 sequences with CD-HIT-EST/PSI-CD-HIT or many other options
Take all the sequences that form one individual cluster and build a multiple sequence alignment (MSA) with MAFFT/ClustalOmega or similar
Take the MSA and build a phylogenetic tree with a Maximum-Likelihood approach like iqtree or similar
Visualize the tree file with Jalview or similar

Of course this is rather general and depending on exactly what you're doing you may want a different workflow and/or different tools. You should also explore the parameter space, do not assume the defaults are necessarily good choices


RDAVIDWebService (R package) manual installation, or alternative tool

That package was deprecated in Bioconductor version 3.13 after several attempts to contact the package author/maintainer.
Your options include reverting back to Bioconductor version 3.12 or using one of the other options for GO/pathway enrichment analyses, e.g. topGO, enrichR, or clusterProfiler.