Normalize RNA seq data from multiple runs for expression analysis

There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization ($(x - \overline{x}) / \sigma$) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that here.
To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.


How can the cell line contribution be estimated from RNASeq data?

There are couple of computational methods which try to do this (I never used them, so no experience):

CellMix, based on sets of marker gene lists
Subset Prediction from Enrichment Correlation, which is based on correlations with subset-specific genes across a set of samples.
Cell type enrichment, which uses our highly expressed, cell specific gene database
Cell type-specific significance analysis using differential gene expression for each cell type

You might have to get some reference expression levels from public databases or papers for some of the methods.

One thing to keep in mind: you cannot really compute the cells proportion, only RNA proportion. If you have a good reason to assume that RNA quantity
per cell is very similar, this is a good proxy for the cells proportion in a tissue.


Proteins don't have ORFs, DNA sequences do.  What DNA sequences do you have?  If cDNA (i.e. transcript/mRNA) then an ORF finder will give multiple ORFs and the coordinates on the full-length input DNA which indicate it's length.


plotMA function issue: can not circle genes

It's likely that the fold-change of the topGene is greater than 5. Consequently, the circle is being drawn and the text produced, but they're outside the bounds of the plot (the bounds don't update if you add a circle and text label). Look at the log2FoldChange column of resLFC1[topGene, ] and increase the bounds to encompass that.


Same transcript coordinates in gtf file, different transcript ID

As Ian explained, these are different transcripts which happen to have the same start and end positions. You have no information on their exonic structure in that file. However, if you look them up at EnsEMBL, you will see:

Transcript Y74C9A.2b.1 :

1   Y74C9A.2b.1.e1  11,499  11,561  -   -   63
    Intron 1-2      11,562  11,617          56
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

This has 4 exons, two of which are non-coding (UTR). The second starts at

Transcript Y74C9A.2b.4

1   Y74C9A.2b.4.e1  11,499  11,557  -   -   59
    Intron 1-2      11,558  11,617          60
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

Granted, the two transcripts are very simlar, they both have 4 exons, the first 2 of which are non coding (UTR). However, the 1st exon of transcript
Y74C9A.2b.1 is 63 nt long and ends at 11,561 while the 1st exon of transcript Y74C9A.2b.4 is 59 nts and ends at 11,557.

So there you go: two (slightly) different transcripts.


Unexpected error using "stdin" in bedtools intersect as a piped command

With v2.27.1 (or any version I ever used) instead of stdin use -
(...) | bedtools intersect -a - -b blacklist.bed -v > out.bam

Be aware that the output of this is a BAM, not SAM file. It would be easier though to simply make a list of regions you want to keep, e.g. the complement of the backlist file with the entire genome (bedtools complement), and then use samtools view -L option to only keep overlapping reads. That saves you the run through bedtools.


What is outliers and how to determine outliers in the following phylogenetic tree?

Distance- and ML-based algorithms using reversible models can't find the root of trees. A classical method to root a tree is to use an outgroup (not outlier), which is a species/sequence known to directly descend from the root of the tree. In your case, it is relatively easy: add a chicken or fish ortholog to your dataset and put the root on the chicken/fish branch.

Outgroup works well for single-copy genes used in species tree construction (e.g. ADH). However, it doesn't always work. The culprit is the "known" part. Say a gene has two copies A and B in vertebrates. Copy A was lost in rodent and copy B lost in primate. If you choose a chicken A gene as the outgroup, the correct tree should be ((primary-A,chicken-A),rodent-B). Without knowing the true history, you may forcefully put the root at the chicken-A branch and build a wrong tree ((primate-A,rodent-B),chicken-A).

There are a few other tree rooting methods. An easy approach is to put the root at the longest branch in the tree, assuming the presence of molecular clock. When the species tree is known, you can root a gene tree by minimizing the number gene duplication/loss events in the history. I generally prefer the latter approach when the relevant information is available.


How to convert VCF format file into Fasta(.fa) with the Reference data

I used FastaAlternateReferenceMaker previously. You will need to download GATK first. Next you prepare your reference genome (reference.fasta) and your vcf file (input.vcf). You will call the GATK function like this:

 java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -V input.vcf \

output.fasta will contain the new fasta with snps inserted at sites specific by the vcf file


STAR-long parameters for aligning RNA ONT reads to genome

I've had great results using minimap2, particularly when combined with a pre-treatment of Canu for error correction (using minimap2 for the read-to-read mapping):

# correct reads
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
   genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
   -p 4T1_BC06 -d 4T1_BC06 \
   -nanopore-raw workspace/pass/barcode06/fastq_runid_*.fastq
# align reads
~/install/minimap2/minimap2 -p 10 \
   -a ~/db/fasta/mmus/ucsc/mmus_ucsc_all_cdna.idx \
   -x splice <(pv all/4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  samtools sort > 4T1_BC06_all_vs_mmusAll.bam

Update: the most recent nanopore basecaller, combined with the most recent minimap2, seems to now do a decent job with mapping, and pre-correction of reads no longer seems necessary. More recently I've been using LAST to map to the transcriptome, and minimap2 to map to the genome. Minimap2 has the ability to use a homopolymer-compressed genome index, which means that the most common consistent error for nanopore reads (i.e. misjudgement of homopolymer length) will not influence the mapping rate.


PAML branch-site model: Identical lnL value between null and alternative models?

The issue with PAML is that it doesn't throw many exceptions (programming jargon for trapping and reporting errors) if something is wrong, generally the "exceptions" are identified in the failure to report an expected value in the output.

Possible explanation fix_omega is not turning the omega parameter on or off it is simply defining what the value of omega will be. If the default value, which is set at .4, but can be changed, is very close to the estimated value the likelihoods will be similar - or even identical.

Goto debugging As computer/data scientists, the very first thing we should do is just check there is run time variation between the options you are assessing via,

time codeml ./file

The thing to remember is fix_omega=0 is where Omega is being estimated from the data, whilst =1 means it will be supplied with a value. There should be longer run time on the "0" option.

Nitty gritty debugging

1 Please check you have

seqtype = 3

If not, then that is your answer, omega was never being set.

2 For debugging only make sure you have,

model = 1

This is global dN/dS for all the data not for a specific clade within the tree, the reason for this is given in the final step below

3, Then try a different NSsites model,

NSsites = 8

This should produce real differences. The NSsites options are,

  0:one w;1:neutral;2:selection; 3:discrete;4:freqs; * 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
  10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1; * 13:3normal>0

w here is for omega

4 The proof of your system is that when you

fix_omega = 1
omega = .4 * try lots of different values, e.g. 0.3, 0.4, 0,7, 1, 1.5, 2

This should produce likelihood variation

Finally

model = 2

Then go to your tree file and make sure that the nomenclature used follows the manual, I think from memory the # is deployed. The output file should change to specify that you have two values for dN/dS from your treefile, thus the output from model=2 should vary from the output of model=1 if not the delineation has been ignored. If you do see two rates, great and now you are into parameter exploration. Parameter exploration is basically the options 1 to 4 described above, plus the stuff you mentioned.

Oh you mentioned using Newick format. I haven't used PAML for a while (dNdS stuff is ancient) but it it only used phylip format when I used it.  I suspect if it got the wrong tree format PAML will throw an exception, but it's worth checking.


How to extract metadata from NCBI's experiment?

You'll want something like:
esearch -db sra -query SRX1596422 | efetch -format runinfo

This will produce a CSV output to the screen with columns containing the meta information available in SRA:
Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR3184375,2016-02-24 00:00:00,2016-02-23 23:14:45,556018223,111203644600,556018223,200,57997,GCA_000001635.1,https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos2/sra-pub-run-7/SRR3184375/SRR3184375.1,SRX1596422,,WGS,size fractionation,GENOMIC,PAIRED,0,0,ILLUMINA,Illumina HiSeq 2000,SRP070623,PRJNA312530,,312530,SRS1306521,SAMN04492320,simple,10090,Mus musculus,MCNT_E1,,,,,female,,no,,,,,THE SCRIPPS RESEARCH INSTITUTE,SRA356784,,public,52A6F836C13F2D15755EDE3E52A18D78,240E7130CFE36405818DA5EE2A10E24C

The more relevant fields are: Subject_ID, Sex, Disease, Tumor, Affection_Status, Analyte_Type, Histological_Type, Body_Site. Note that not all of these are available for that particular sample.


Is it possible to create a DESeqDataSet with a user-provided design matrix?

Provide rank sufficient design to DESeqDataSetFromMatrix and then use your custom model matrix in DESeq. In essence:

dds = DESeqDataSetFromMatric(counts, s2c, design=~batch)
design <- model.matrix(~strain+batch, s2c)
design = design[, -9]
DESeq(dds, full=design)

See this thread on the bioconductor site for details.


Let $n$ be the number of rows (in your case $n=5$) and $|O_j|$ the cardinality of $O_i$.

We can calculate the number of compatible columns as sum of the number of $O_i \subseteq O_j$, $O_j \subseteq O_i$ and $O_i, O_j$ disjoint.

number of $O_i \subseteq O_j$:
Here we need to find the number of supersets of our $O_i$. That means we need to find the number of sets which contain the indices in $O_i$ and for each index not in $O_i$, we have two choices: It is contained, or it is not contained in $O_j$. Therefore the number of supersets is $2^{n-|O_i|}$, which is 2 to the power of elements not in our set.

number of $O_i, O_j$ disjoint:
Here all elements of $O_i$ are fixed not to be in $O_j$ and for each index not in $O_i$, it can again be in $O_j$ or not. That means that the number of possible sets is again $2^{n-|O_i|}$.

number of $O_j \subseteq O_i$:
Here we need to find the number of subsets of $O_i$ and this is the number of elements in its power set, which is $2^{|O_i|}$.

So finally, this means that the number of compatible columns can be calculated as: $2^{n-|O_i|} + 2^{n-|O_i|} + 2^{|O_i|} - 2$.

Why -2? We have counted 2 sets twice:

The empty set, which is a subset of $O_i$ and also counted in the number of disjoint sets
and $O_i$, which we counted in $O_i \subseteq O_j$ and in $O_j \subseteq O_i$

This means for your examle, we have $16 + 2 + 2 - 2 = 18$ compatible sets.

Note: The all zero vector (which corresponds to $O_i=\{\}$) is a special case for which the formula above does not work, as all disjunct sets and all supersets are counted twice. For this special case, the number of compatible columns is simply $2^n$.


How to seperate the genus into different group based on th ANI calculation?

ANI = average nucleotide identity.

This will give you a 2x2 matrix. You reconstruct a phylogeny using clustering, the most important being neighbor-joining, UPGMA is very occassionally used.

The statistic is clustering algorithm and its done by pairwise reconstruction of sister groups. Once a sister-group is defined (usually on nearest-neighor both in neighbor-joining and UPGMA) the genetic divergence from the hypothetical node is established either using the mean (UPGMA) or additive criteria (neighbor-joining). Okay enough theory ...

You get Mega X and find out the import format for a 2x2 distance matric. Mega X however will calculate this for you. You then go to the build tree menu heading and click "neighbor-joining" (or "distance"). The tree will very quickly automatically appear.


I have used TER to break the long bonds of a chain in my PDB

Normally you have three options in a PDB:

Add a TER record
Add the  OXT atom, if your chain was cleaved or cleaved spontaneously (i.e. Asp-Pro sequence), the C-terminal would sill have a terminal oxygen. i.e. making it a terminal residue. If you have protons, you'll need to add also 3H .
Assign a different chain

If you are using forcefields as is your case, the second is the only choice. The C  atom and the N  atoms have a type that does not match because they have the wrong number of bonds —hence the type error.


Where to download JASPAR TFBS motif bed file?

You could make your own such file with FIMO, the JASPAR MEME files available via the MEME download site, and your genome build of interest (per-chromosome FASTA files from UCSC goldenpath, say), e.g., for assembly hg19 and UCSC chromosome naming scheme:
$ for chr in `seq 1 22` X Y; do echo ${chr}; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/chromosomes/chr${chr}.fa.gz | gunzip -c - > hg19.chr${chr}.fa; done

For each chromosome (omitting the for loop code for brevity), you can use the fasta-get-markov in the MEME suite to generate a 5th-order Markov model of kmer frequencies:
$ fasta-get-markov -m 5 hg19.chr${chr}.fa > hg19.chr${chr}.5th_order_background_model.txt

Here is one way to download a non-redundant, vertebrate JASPAR database (MEME-formatted):
$ wget http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.txt

To run a FIMO query with the 5th-order background model and JASPAR database, using a p-value threshold of 1e-4:
$ fimo --verbosity 1 --bgfile hg19.chr${chr}.5th_order_background_model.txt --thresh 1e-4 --text JASPAR2018_CORE_vertebrates_non-redundant_pfms_meme.txt hg19.chr${chr}.fa \
  | tail -n+2 \
  | awk -vOFS="\t" -vFS="\t" '{ print $3, $4, $5, $1":"$2">"$10, $8, $6 }' \
  | sort-bed - \
  | starch --omit-signature - \
  > hg19.chr${chr}.bed.starch

The output of a FIMO search is a text file which can be modified into a sorted BED file via sort-bed, an example of which is shown above.
To get a faster answer, each per-chromosome query can be put onto its own node on a computational cluster.
Using a p-value threshold of 1e-4 will result in very large files (when uncompressed). Using tools to compress BED files is useful here. Use of BEDOPS starch to compress BED files (as shown above, for instance) will give better compression performance than bgzip with per-chromosome random access and integration with BEDOPS bedops and bedmap CLI tools. The bgzip option offers integration with tabix and offers finer-grained random access, but it will use more disk space and doesn't integrate with BEDOPS tools.


Generate SNP/indel annotation in Arabidopsis

The Ensembl Variant Effect Predictor will annotate your A. thaliana variants:
http://plants.ensembl.org/Oryza_sativa/Tools/VEP

The web based tool will use the reference genome assembly and gene annotation used in the Ensembl browser:
http://plants.ensembl.org/Arabidopsis_thaliana/Info/Annotation/#assembly

However, if you use the command-line version of VEP, you can provide your own reference genome assembly and gene annotation files to be used as a reference.

You can find all of the documentation here:
http://www.ensembl.org/info/index.html

Best wishes

Ben
Ensembl Helpdesk


Ancestry of the coronavirus 2019-nCov, WuHan city, China

Overview
The central focus of the tree is to highlight the key biological concern of the new coronavirus, 2019-nCov. The key concern is the genetic similarities to SARS epidemic, and relates to the SARS receptor.

SARS background SARS is endemic in bats (your BioRxiv tree partly shows that and this tree definitely shows it) and in the 2002 epidemic infected civet cats which then infected humans. More importantly it is then transmitted from one human to another. This is a consequence of the SARS receptor been able to exploit the cellular receptors in the respiratory tract of bats, civet cats and humans and use this tissue as a site of replication. The concern is exacerbated because the tree here shows SARS independently infected humans on two separate occasions, suggesting the cross-over has an underlying genetic basis, which is part of the "SARS lineage".

BioRxiv 2019-nCov tree Your tree shows the 2019-nCov has a recent common ancestor with SARS in comparison to the rest of the betacoronaviruses. This therefore provides some circumstantial evidence the receptor mechanism and ability to frequently crossover from bats and ultimately transmit human to human could be shared with SARS. At an evolutionary level (which is what ANY tree is supposed to reconstruct) it raises the question which common ancestor of SARS did the "SARS receptor" originate. I'd need to draw a diagram to better demonstrate this point, but I hope you get the idea.

The "influenza receptor" The analogy is with influenza virus and the entire epidemiology of influenza hinges around is sialic acid receptor of influenza and its ability to bind on to the cellular receptors in the upper/lower respiratory tract of birds, pigs and humans. The hypothesised mechanism is called the mixing vessel theory and is the classic epidemiological understanding of how new pathogenic influenza pandemics occur. If you replace "birds and pigs" with e.g. "bats, civet cats and humans" you get the idea why this could be scary. We don't the intemediate host of 2019-nCov ... but I speculate there must be one, unless eating Cov-infected bats is common in China.

Technical details of the BioRxiv tree The tree is a nicely diverse selection of the beta-coronaviruses. The authors have rooted the tree using the outgroups of delacoronaviruses and gammacoronaviruses, so it is a good robust selection of outgroups which can be used to correctly identify the direct of evolution of betacoronavirus divergence. In tree theory (its formal name is phylogenetics theory), extensive rooting is good and minimises artefacts.

One of the earliest members of the betacoronaviruses to diverge is MERS (Middle Eastern respiratory syndrome), which form a single "clade" (all viruses share a unique common ancestor) and these represent around 50% amino acid divergence from SARS/2019-nCov. The selection of MERS in the tree from both camels and humans looks good. The selection of other betacoronaviruses looks great, I wasn't aware of the "ruminant clade"  at all involving buffalo, cow etc .. infections and there was an associated human infection. There are loads of bat isolates of the betacoronaviruses throughout the tree, but we eventually arrive at the "SARS clade". The authors show that 2019-nCov is an outgroup to the SARS clade and shows a close relationship to one (BioRxiv tree) or two (this tree) bat isolates. Looking at the precise SARS clade (better shown in this tree) we can see loads of bat virus associating with SARS lineages. We therefore assume the reservoir to SARS and likely 2019-nCov is bat, moreover that the single ancestor to both viruses was a bat (it is called a parsimonious hypothesis).

The one thing the BioRXiv tree omits, for example with regards this tree, is the diversity of SARS and in particularly the two independent origins of SARS, which is a weakness of their analysis, particular if this fed into downstream analysis. It is not to say the authors were wrong, but it was uncool.

Word of caution my understanding is the divergence between SARS and 2010-nCov is around 15% and this is a quite a large amount of genetic divergence, even if both viruses share a most recent common-ancestor, form part of the same clade and receptor peptide motifs. Nevertheless there is sufficient amino acid divergence to generate notable differences in epidemiology, clinical symptom and transmission.

MERS, SARS and 2019-nCov mortality rates It is worth noting that within the betacoronaviruses the mortality rates between difference clinically important viruses is very different. MERS has a mortality rate of 40-60%, SARS is around 10% but 2019-nCov is <2.5%. The 2019-nCov mortality rate is still important given it has infected more and could potentially infect alot more people than SARS.


Counting the number of paralogues for mouse genes gives me the wrong frequency in R

I'd use dplyr for this.

First, mouse paralogues contains no values for mmusculus_paralog_associated_gene_name that are NA:

library(dplyr)
mouse_paralogues %>%
  filter(is.na(mmusculus_paralog_associated_gene_name)) %>%
  nrow()

[1] 0

But there are values which are empty strings, so you can filter for those if you wish.

count_mouse_paralogues <- mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name != "") %>%
  count(external_gene_name)

Confirm that this gives the expected result (n = 4) for Plin4:

count_mouse_paralogues %>%
  filter(external_gene_name == "Plin4")

# A tibble: 1 x 2
  hsapiens_homolog_associated_gene_name     n
  <chr>                                 <int>
1 Plin4                                     4

And check mouse_paralogues:

mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name == "PLIN4")

  hsapiens_homolog_associated_gene_name external_gene_name mmusculus_paralog_associated_gene_name
1                                 PLIN4              Plin4                                  Plin3
2                                 PLIN4              Plin4                                  Plin5
3                                 PLIN4              Plin4                                  Plin2
4                                 PLIN4              Plin4                                  Plin1


I couldn't find a package neither in Github nor in Bioconductor. Neither could I find any paper about a package linking to bioDBnet. So I don't believe there is any (R) package at the moment which links to biodbnet in order to retrieve the IDS.
However it has a couple of APIs so you could build one or use them to retrieve the information you want.


Multiple sequence pairwise global alignment against reference sequence

Q1

From this phrase, "The end goal is to identify every unique sequence", I am not 100% sure of what you mean by "every unique sequence". It sounds like maybe you want to group alignments of homologs and visualize them together?

If the proteins are similar enough, you might also benefit from trying BLAST to look at similar regions.

Otherwise MAFFT, Clustal, or MUSCLE are all appropriate tools to generate the alignments that you need. However I also do not know of any software aside from Geneious that has the option to align many proteins to a reference. If you were working with DNA I'm sure the stock suggestion would be to use a read aligner like bwa mem or minimap2 and to use a bam file.

Since I don't know of any software to do exactly what you would like, I wrote a small python script to run MAFFT many times for each sequence in a fasta after the first sequence, the reference.

Q2 - MAFFT/Python

In any case, here is a bit of python code that performs a MAFFT alignment of all of the sequences in a fasta file against the first sequence in the fasta file, which you have designated as the reference.

Here is the fasta file from Wikipedia. I called it seqs.fasta:

>reference
MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*

>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
DIDGDGQVNYEEFVQMMTAK*

>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY

Here is a script, mafft_ref.py. Run it with python mafft_ref.py seqs.fasta:

from Bio import SeqIO
from Bio.Align.Applications import MafftCommandline
import tempfile

fasta = sys.argv[1]

reference = 0
count = 0
with tempfile.NamedTemporaryFile() as temp:
    for record in SeqIO.parse(fasta, "fasta"):
        if count == 0:
            reference = record
        else:
            SeqIO.write([reference, record], temp.name, "fasta")
            mafft_cline = MafftCommandline(input=temp.name)
            stdout,stderr=mafft_cline()
            print(stdout)
        count += 1

Here is the output of python mafft_ref.py seqs.fasta printed to the terminal. Use python mafft.py prots.pep > output.txt to save the output to output.txt.

>reference
------------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQ
VSLRDLFDRA--------------VMVSHYIHDLSS-----EMFNEFDKRYAQGKGFITM
ALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHLVTEVRGMKGAPDAIL
SRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDEDARYSAFYNLLH
CLRRDSSKIDTYLKLLNCRIIYNNNC
>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKEL------------GTVMRSL---GQNPTE
AELQDMINEVDADGNGTIDFPEFLTMMARKMKDTDSEEEIREAFRVFDK---DGNGYISA
A---------------------------------------ELRHVMTN-----------L
GEKLTDEEVD-------EMIREADIDG-------------------DGQVNYEEFVQMM-
-----TAK------------------

>reference
------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVS----
--LRDLF--------------------DRAVM----VSHYIHDLSSEMFNEFDKRYAQGK
GFITMALNSCHTSSL-------PTPEDKEQAQQTHHE-------------VLMSLILGLL
RS--WNDPLYHLVTEVRGMKGAPDAILSRAIEIEEE--------------NKRLLEGMEM
IFGQVI-------------------P-------------------GAKETE-PYPVWSGL
PSLQTKDEDARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC
>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGY----VLPWG----QMSFWGA
TVITNLFSAIPYIGTNLVEWIWGGFSVDKATLNRFFAFHFILP-----------------
-FTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLGLLILILLLLLLALL
SPDMLGDPDNHM----------PADPLNTPLHIKPEWYFLFAYAILRSVPNK--LGGVLA
LFLSIVILGLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQM
ASI----------LYFSIILAFLPIAGXIENY--------------

References:
Some credit of mafft_ref.py to a biostars answer found here.


Why is FPKM still used for gene expression studies?

Is it ever meaningful to use FPKM values for analyzing across samples?

One should never use FPKMs for anything important. They can occasionally be useful for plotting, but even in that case one needs to construct the FPKMs from properly normalized data and not use the original definition.

It's becoming ever rarer to see people using FPKMs these days, thankfully.


how to easily make complex Euler/Venn diagrams in R?

Since you are in R, consider using the UpsetR library.  It doesn't make Venn diagrams but it helps to visualize overlaps between any number of groups.
http://gehlenborglab.org/research/projects/upsetr/#:~:text=UpSetR%20is%20an%20R%20package,based%20on%20groupings%20and%20queries.


Phyre2 vs ITasser, completely different models generated

I'm less familiar with Phyre, but I-TASSER is a really sophisticated system that takes the results of a search using multiple threaders and plugs them into an ab initio simulation which tries to minimize the energy of the models by sampling many possible 3D conformations, which I don't think Phyre does.

https://en.wikipedia.org/wiki/I-TASSER#/media/File:I-TASSER-pipeline.jpg

Compare with a similar workflow schematic for Phyre:

http://www.nature.com/nprot/journal/v10/n6/images/nprot.2015.053-F1.jpg
http://www.nature.com/nprot/journal/v10/n6/fig_tab/nprot.2015.053_F2.html

Structure prediction still has a long way to go, and you'll always get better results if there are close homologues available in the PDB, but given the consistent high performance of I-TASSER in CASP I would treat those results as more significant. That said, it can't hurt to consider multiple answers.

Edit: included blmoore's link to second half of Phyre protocol schematic


Canu assembly not making a single consensus?

Looking that the link you provided only 30% of the reads appear to be e. coli.  I would suggest filtering the reads for those that align to e. coli and see where that gets you.  Did the original publication describe the method they used for data cleaning?

A:

https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR3385327

Unidentified reads: 68.82%
Identified reads: 31.18%
cellular organisms: 30%
Bacteria: 30%
Proteobacteria: 28.09%
Gammaproteobacteria: 19.7%
Enterobacterales: 13.3%
Enterobacteriaceae: 12.54%
Vibrionales: 0.01%
Viruses: 1.17%

B:

Without a good assembly canu won't be able to determine if the assembled genome is circular.  I would ignore this until you have a good example.


installing multiple bioconductor packages at once

I always do something like this:
bio_pkgs <- c("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData",
          "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

# install:
BiocManager::install(bio_pkgs)

# load all at once
invisible(lapply(pkgs, function(x) library(x, character.only=TRUE)))


How can I apply proportional (p) distances (Nucleotide) using bioPython

Based on the sites you mention, it seems that p-distances are pairwise distances between aligned sequences obtained by counting the number of differences and dividing by the length of the alignment. This is the proportion of differences between two sequences.

I'm not sure exactly what you mean by "add p-distance model in my code", but here is how you could compute those distances from an existing alignment (here I use as an example an alignment I have in fasta format):

from Bio import AlignIO

def p_dist(seq1, seq2):
    ali_len = len(seq1)
    assert ali_len == len(seq2)
    return sum(0 if nuc1 == nuc2 else 1 for (nuc1, nuc2) in zip(seq1, seq2)) / ali_len

ali = AlignIO.read("most_targeted_histones/transcript_alignments/H4_top4_loss_spliced_NT.fa", "fasta")
nb_seq = len(ali)
distances = {}
for i in range(nb_seq - 1):
    rec1 = ali[i]
    j = i + 1
    while j < nb_seq:
        rec2 = ali[j]
        distances[(rec1.id, rec2.id)] = p_dist(rec1.seq, rec2.seq)
        j += 1

print(*distances.items(), sep="\n")

Output:

(('ZK131.4_spliced', 'ZK131.8_spliced'), 0.0)
(('ZK131.4_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.4_spliced', 'ZK131.1_spliced'), 0.0)
(('ZK131.8_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.8_spliced', 'ZK131.1_spliced'), 0.0)
(('T10C6.14_spliced', 'ZK131.1_spliced'), 0.10897435897435898)

If you want to build a tree from those distances, I'm not sure you can do this with BioPython. It seems however that there are some possibilities with DendroPy.

You would need to write the distances in a suitable format, and to choose the method between UPGMA and neighbour-joining.


Downloading dataset from SRA (SOLiD Platform)

You have to use the accession name without the .sra suffix.

abi-dump -A SRR1657115


This is almost right, but it's missing a little context. In ChIP-seq, it's typical to shift the positive and negative reads because otherwise the peaks end up offset from each other. The amount of shift is determined by trying different shifts and measuring how well the read densities correlate between positive and negative strands. In the question, the first sense talks about a peak on the plot of read density vs genomic position. The second has different x and y axes because it's talking about the cross-correlation between strands. X is the offset between positive and negative strands, and y is the correlation observed at that offset. Look at figures 4D,E in this ENCODE paper for visuals.
Landt, S. G., Marinov, G. K., Kundaje, A., Kheradpour, P., Pauli, F., Batzoglou, S., ... & Chen, Y. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome research, 22(9), 1813-1831.
https://genome.cshlp.org/content/22/9/1813.full.pdf


Convert BAM to properly paired FASTQ files

As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of samtools collate.

… | samtools fastq -F 0x900 -@ 48 \
    -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -

Your samtools fastq command is not doing anything to siphon off these singleton reads. You should add -s /dev/null to get rid of them.

See the discussion in samtools/samtools#874 (especially the part starting here). I have not looked at how the recent changes from jkbonfield, which will land in an upcoming samtools version (after 1.9), might affect this.


Extract autodocked protein-ligand connections programatically

With pymol as a python module (conda install -c schrodinger pymol, not the GUI one), it is very easy.

import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load('protein.pdbqt', 'protein') # the second is the name within pymol instance.
     pymol.cmd.load('ligand.sdf', 'ligand')
     pymol.cmd.save('neighbours.pdb','ligand expand 3') # get atoms 3 Å around

Done!
However, if you are doing stuff in python in the first place you might want the last line as a string: pymol.cmd.get_pdbstr('ligand expand 3'). Changing the selection 'ligand expand 3' to 'byres ligand expand 3' will save the full residues.

Alternatively you can get the atom coordinates themselves etc. with it. To see what attributes Atom has you have you can do this:

for atom in pymol.cmd.get_model('ligand expand 3').atom:
    print(dir(atom))
    break

Which will give a large amount of attributes, such as atom.coord, which is a 3-element list.


Does GISTIC (v 2.0) estimate amplified/deleted probabilities on a single sample basis?

There doesn't seem to be much differences between GISTIC 1.0 and 2.0 as it says:

  As with GISTIC 1.0, we obtain P-values for each marker by comparing the score at each locus to a background score distribution generated by random permutation of the marker locations in each sample

But on the supplementary material of the GISTIC 1.0 there is a more detailed explanation of the method. See the "Stage 2" section:

  Second, we compare these G-scores to the distribution of scores expected if only random aberrations were observed. This distribution can be determined by rescoring the genome after permuting marker locations within each sample; we instead derive a semiexact estimate.

Furthermore in another section ("Stage 2: Aggregation of Data from Different Tumors to Differentiate Between Driver and Passenger Aberrations") it says:

  To determine which of the aberrations identified in Stage 1 are likely to represent driver events, we aggregate the data from all tumors used in the analysis to generate summary scores for amplifications, deletions, and LOH. The statistical significance of each score is determined by comparison to the distribution of scores obtained by all permutations of the data (using a semiexact approximation), with correction for multiple hypothesis testing.

However, the relevant section (of the supplementary materials) seems to be "Null Hypothesis Generation: An Analytic Derivation of the Null Distribution", where it describes the semiexact approximation used.

In a supplementary file it describes it as "Generate all permutations of SNP labels within each sample to simulate datasets with random aberrations"

Conclusion:

It seems that the background probabilities are calculated within sample.

I can't say if it matters how the sets are constituted but I would say that the broader the sets are, the better estimation of the structural variations it will perform but the same estimation of the background will be done.

Ultimately you can check the code of the program, or test with 10 samples and replace one of them to see if the results are changing accordingly.


Maybe take a look at workflow management systems, like

snakemake
Nextflow
Toil
Cromwell

Admitted, each of those adds their additional learning curve. However, especially snakemake is not so much more complicated than a standard bash workflow. Those systems do help a lot organize and orchestrate your workflows


Outline of solution:

get this file: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt
filter for lines having "Complete Genome" in column 12
filter for lines having a taxid (column 6) that is a descendant of taxon id 1117 (phylum Cyanobacteria), you can use nodes.dmp from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
The file download URL is the concatenation of column 20 and the last field of column 20, after separating it at '/', followed by appending '_genomic.gbff.gz'.  For example: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/007/365/GCF_000007365.1_ASM736v1/GCF_000007365.1_ASM736v1_genomic.gbff.gz

Edit:
Since I will need something like that soon, I made a Perl script for downloading genomes.
For Cyanobacteria, you would do ./download-refseq-genomes.pl 1117


How to calculate various properties of a protein structure per atom using PDB2PQR and Rosetta tools (or any other tools)?

Per Atom
You are correct. Rosetta scorefunction does not store any per atom data. The scoring operates at the per residue level. Whereas each atom has its coordinates and properties in full atom mode, in Pyrosetta it is clear that an atom is a just part of a residue and every operation is applied at the residue level. It's a team effort: the functional group of the side chain works as one and in most cases two atoms have the same role (e.g. OE1 and OE2 in glutamate behave identically).
So, the real question is whether you really care about the atom level in a classical mechanical setting?
Formal charges
A wee parenthesis, we need to make a difference between formal charge and Marsilli-Gasteiger partial charge —one is an integer and the other is a float. Both require a stated protonation state in most applications —say aspartate vs. aspartic acid.
PDB2PQR formats to a PQR file for the ligand docking application autodock. It does not assign how a residue is protonated.
In the case of Rosetta, the partial charge depends on the 4 column of the topology (params file) of the amino acid, i.e. it is fixed for all amino acids —Rosetta does not implement a Drude particle to make something polarisable. I very strongly suspect this is the same in PDB2PQR.
Therefore, something like the Advance Poisson-Boltzmann Solver may give more interesting data in terms of partial charges, but also requires the uses to specify the protonation state.
A small parenthesis is histidine, which can adopt four forms. Two tautomers at neutral pH: single proton on either the epsilon (HIE in Amber) or delta nitrogen (HID), while it can be deprotonated (anion) or diprotonated (HID, cation). Rosetta switches between these two automatically.
Properties
You list a few properties of interest, which all seem to be based on the residue a priori of structure.
For example:

hydrophobicity, unless you have a polarisable model (Drude particle) this is structure independent. RSA, relative surface area, is a must better indication (per residue), which then can get classified as core or surface based on this. Per atom would require just a wee tinker.
aromaticity, this is also structure independent, but π-π stacking or π-cation interactions are noted at the residue level.
hydrogen bond acceptor/donor. Whether an atom is a donor or acceptor is actually defined in the atom type used in the topology file. Whether that atom fulfils that bonding is a different matter. If you really really wanted to get per atom values of this, in Rosetta you could exclude the backbone atoms you could mutate to alanine without minimising and see if it's maintained. So differentiating between backbone and sidechain is fine, but per atom is a different matter.

A structure-dependent property is how flexible the residue is, as measured as an RMS. This can be calculated per atom. This is akin to B-factors, which is how big is the "blur" of the atom. However, if you get the NZ atom on a lysine that is switching between two or more positions you get a partial occupancy. If this atom is basically hand-waving along an axis in an arc you get an anisotropic B-factor. In other words it does not vibrate in all 3 cartesian axes equally. So a B-factor is a terrible approximation. An RSMD of an ensemble is even worse.
In Rosetta, you can a quick vibration using Backrub, which is very quick at telling you how much motion does an atom experience.
Here be dragons

Warning: As mentioned, it takes two atoms to resulting in a force, hence why per atom scores are a bad idea.

However if for some exceptional reason this is required, say
figuring out which group on a ligand is causing a clash, in pyrosetta it is possible:
import pandas as pd
import pyrosetta
pyrosetta.init()
pose = pyrosetta.pose_from_sequence("ELVISISALIVE")
target_resi = 9
residue = pose.residue(target_resi)
scores = []
for i in range(1, residue.natoms() + 1):
    iname = residue.atom_name(i)
    for r in range(1, pose.total_residue() + 1):
        other = pose.residue(r)
        for o in range(1, other.natoms() + 1):
            oname = other.atom_name(o)
            score = pyrosetta.toolbox.atom_pair_energy.etable_atom_pair_energies(residue,
                                                                                i,
                                                                                other,
                                                                                o,
                                                                                scorefxn)
        scores.append(dict(zip(['residue_index_A', 'atom_name_A', 'residue_index_A', 'atom_name_B', 'lj_atr', 'lj_rep', 'fa_solv', 'fa_elec'],
                         [target_resi, residue.atom_name(i), r, other.atom_name(o), *score])))

df = pd.DataFrame(scores)
df.pivot_table(index='atom_name_A',
               aggfunc='sum')[['lj_atr','lj_rep','fa_elec','fa_solv']]


Python module for fetching NCBI id for a list of species

Biopython has an efetch frontend.

from Bio import Entrez
Entrez.email = "michaelg@example.ac.uk"
handle = Entrez.efetch(db="nucleotide", id="KR052012.1", rettype="gb", retmode="text")
print(handle.read())

OUT

LOCUS       KR052012               10714 bp    RNA     linear   VRL 09-MAY-2015
DEFINITION  Dengue virus strain Hb33/CHN/2014, complete genome.
ACCESSION   KR052012
VERSION     KR052012.1
KEYWORDS
.... etc

You will likely need to perform an EFetch search and then extract out the accession numbers for import (as done above)