Basecalling process carried out by full_1dsq_basecaller.py/.exe

My understanding is that the 1D$^2$ process is carried out in three stages:

Carry out a 1D calling of reads (i.e. what the standard Albacore does at the moment)
Link reads that are similar to each other (in reverse-complement orientation), from the same channel, and have a very short pore transition time
Create a 2D consensus base call, using information from the two linked reads

The two outputs from the basecaller are the 1D calls from step 1, and the consensus calls from step 3. Ideally just the 2D consensus calls should be used for assembly (assuming it represents a random distribution of reads), but  if that gives too low coverage (or is not random enough), then the non-converted reads from the 1D calls could be included as well.

Because the linking is not a physical link, it's not a good idea to use 1D$^2$ for amplicon reads; the likelihood of two complementary reads from different templates appearing one after another is too high.


Alternative splicing positions on protein

If you are considering the human genome as "data" (not mRNA), I would run your single peptide via tblastn which will take your protein(s) translate them to nucleotides and run a blast against Genbank. The splice sites will be immediately apparent because they force indels in your query. If there are no indels there are no introns (splice sites). If you want to process several proteins, one solution is to use eblast, a command line version of tblastn  (Anaconda package). You can confirm this against a full genome annotation (note, I don't work on humans).

Alot of work has been done on splice sites, a friend worked on this, but I think you just want to focus on a restricted number of target peptides in your mass spec. output.

mRNA should be avoided IMO if you are looking for splice sites, i.e. intron boundaries.


How to sort multiple FASTA files based on their content?

That's quite easy with cmd:
@echo off
for /l %%i in (1,1,8) do (
  for /f "delims=" %%f in ('findstr  /mc:" segment %%i " *.txt') do (
    ECHO move "%%f" "segment %%i\"
  )
)

findstr /mc:"string" *.txt searches each matching file for the string and outputs the filename, if the string was found. The for /f loop around catches those filenames, so you are able to move them. The for /l loop just iterates from 1 to 8.
See for /? and findstr /? for more information about those commands.
You probably want to adapt the filemask *.txt to your needs.
Note: I "disarmed" the move command for troubleshooting. Remove the ECHO when you are sure it does what you want.


Difference between BWA-backtrack and BWA-MEM

TL;DR:

BWA-backtrack is based on backtracking.
This approach is appropriate only when the dissimilarity between the reads and the reference is low,
or when you want to find all best hits or enumerate all possible alignments up to a specified
number of errors.

In all other situations, BWA-MEM is preferable as it can,
thanks to its sophisticated strategy based on maximum exact matches,
deal with errors better and also automatically switch between local and
global alignment modes.

Long description:

I would like to provide some algorithmic insight since I believe that it
can be very useful in this case.
Both BWA-backtrack and BWA-MEM use the same indexing strategy
(heavily relying on BWT-index), but the actual algorithms are quite different.

BWT-index (and also other full-text indexes such as
suffix arrays or suffix trees)
can easily find exact matches (imagine Ctrl+F-like search in a text editor),
but any differences between the read and the reference,
such as sequencing errors or genomic variants,
make the situation complicated. One then needs to
somehow transform inexact matching to exact matching,
and all three BWA mappers (note that there exists also BWA-SW, but it is deprecated)
use quite different strategies.

BWA-backtrack looks for substrings of the reference, which would be similar
to the entire read (end-to-end) using an algorithm
called backtracking.
First, it searches occurences of the read without any "corrections".
If nothing found, it consideres all possible single edits; then two edits, etc.
To make mapping efficient, one usually wants to stop with the first
found alignment as this would be the best one.
When required, it is also possible to find the other equally good alignments
or to enumerate all alignments up to some edit distance or withing some divergance rate
(see the -N option of BWA-backtrack).

It turns out that
the time required for finding an alignment
can be exponential in the number of errors, which is probably the main problem of
backtracking-based approaches.
To prevent huge overheads due to dissimilar reads,
one needs to limit the number of allowed errors to some reasonable number
(see maxDiff in the BWA man page)
and consider the other reads unaligned.
In the case of BWA-backtrack,
the minimum required identity level is ~97% with the default options (see the -n option).

In fact, the algorithm is more complicated and uses various heuristics such as
seed-and-extend or Z-dropoff in order to make the computation fast enough (at the price of lower accuracy). If you are interested
in more details, all these tricks are well described in the paper.

BWA-MEM uses quite a different strategy. It detects long exact matches
between the read and the reference, and then chains them into local or global alignments,
based on what is more appropriate in that specific case.
Such an automatic local-global switching can be very powerful and
BWA-MEM works well with various types of data
(short reads, long reads, low error rates, high error rates, etc.).


Merging sequencing data for ChIP-seq experiments

I don’t think it matters. Both are easy to merge (BAM via samtools merge, and (gzipped) FASTQ via cat), and neither method has specific disadvantages, unless your FASTQ files are sorted for some reason (but they generally shouldn’t be).

One advantage of keeping the FASTQ files separate is that it makes it slightly easier to parallelise the mapping step: just run the mapper in parallel on the separate FASTQ files. Although bowtie has an option (-p) for this, throughput from that is slightly worse than running the mapping on split files.


A de novo assembly with spades should work as well as a mapping with bowtie2.
While performing a de novo assembly you should also be aware of possible contamination which you may need to remove. By using a mapping this may not be necessary.
When you use a mapping, you must be aware that you may miss relevant and diverged regions because of the mapping to a reference strain. So a mappping would be preferrable if your target genome is closely related to your reference genome.
I'd like to give you the advice to use some kind of distribution system like Anaconda which helps you to keep your system clean and make it easier to reproduce your steps.
Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :
De novo assembly
1. Quality Check ( you already did this )
2. Run Spades
spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &> k12.log

This command will start an assembly using 4 threads and 8 GB RAM. You can scale it to whatever your system is capable of or what you'd like.
You now have a folder called k12 which holds the assembly information.
To get some information about the contigs you can use the emboss toolkit.
infoseq k12/contigs.fasta

fasta::contigs.fasta:NODE_859_length_78_cov_1.000000 -              NODE_859_length_78_cov_1.000000 -              N    78     43.59
fasta::contigs.fasta:NODE_860_length_78_cov_1.000000 -              NODE_860_length_78_cov_1.000000 -              N    78     44.87
fasta::contigs.fasta:NODE_861_length_78_cov_1.000000 -              NODE_861_length_78_cov_1.000000 -              N    78     43.59

This is the end of the contigs.fasta. So you see contigs.fasta contains 861 contigs. And you also see thate ther length and coverage is quiet low. So you need to remove this low coverage regions.
3. Remove low quality contigs
There are multiple ways how you could achieve this. I'd like to filter sequence smaller than 500 bp. the tricky part is to identify the coverage threshold. Contigs with a low coverage could be some artifacts or contamination.
4. Check for possible contamination
So to check for possible contamination you can use kraken2. But it needs some disk space. I am using here the minikraken database
kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &> kraken.log

The ecoli_kraken.report shows an overview of how the contigs could be classified to availible species from the database.
Based on this results you can also remove some contigs.
5. Run Prokka
Now after you removed unwanted contigs, you can run prokka to get the annotation.
prokka --cpus 4 --outdir prokka k12/contigs.fasta &> prokka.log

And you will find the information in the prokka folder
Mapping
The mapping process is a little bit easier. You can use bowtie2 and samtools
bowtie2 -x ecoli -1 SRR6940089_1.fastq -2 SRR6940089_2.fastq -p 4 | samtools view -@ 4 -O BAM -o ecoli.bam

158216 reads; of these:
  158216 (100.00%) were paired; of these:
    35696 (22.56%) aligned concordantly 0 times
    120047 (75.88%) aligned concordantly exactly 1 time
    2473 (1.56%) aligned concordantly >1 times
    ----
    35696 pairs aligned concordantly 0 times; of these:
      4458 (12.49%) aligned discordantly 1 time
    ----
    31238 pairs aligned 0 times concordantly or discordantly; of these:
      62476 mates make up the pairs; of these:
        56512 (90.45%) aligned 0 times
        5365 (8.59%) aligned exactly 1 time
        599 (0.96%) aligned >1 times
82.14% overall alignment rate

You now also get some information about the similarity to your reference sequence.
Now you need to work a little bit more on the BAM file. You need to sort and index this.
samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam

Now you can create a consensus sequence using bcftools and seqtk
bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N > ecoli_consensus.fasta

This gives you now a sequence where unknown bases are marked as N.
As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file bio.yaml
name: bio
channels:
  - bioconda
  - defaults
dependencies:
  - aragorn=1.2.38=h470a237_2
  - barrnap=0.9=2
  - bcftools=1.9=h4da6232_0
  - bedtools=2.27.1=he941832_2
  - blast=2.6.0=boost1.60_0
  - bowtie2=2.3.4.3=py36h2d50403_0
  - emboss=6.6.0=h6debe1e_0
  - hmmer=3.2.1=hfc679d8_0
  - infernal=1.1.2=h470a237_1
  - kraken2=2.0.7_beta=pl526h2d50403_0
  - libdeflate=1.0=h470a237_0
  - libidn=7.45.0=2
  - minced=0.3.2=0
  - parallel=20160622=1
  - perl-app-cpanminus=1.7044=pl526_1
  - perl-bioperl=1.6.924=4
  - perl-carp=1.38=pl526_1
  - perl-constant=1.33=pl526_1
  - perl-data-dumper=2.161=pl526_2
  - perl-encode=2.88=pl526_1
  - perl-exporter=5.72=pl526_1
  - perl-extutils-makemaker=7.34=pl526_2
  - perl-file-path=2.15=pl526_0
  - perl-file-temp=0.2304=pl526_2
  - perl-parent=0.236=pl526_1
  - perl-threaded=5.22.0=13
  - perl-xml-namespacesupport=1.12=pl526_0
  - perl-xml-parser=2.44=pl526h3a4f0e9_6
  - perl-xml-sax=1.00=pl526_0
  - perl-xml-sax-base=1.09=pl526_0
  - perl-xml-sax-expat=0.51=pl526_2
  - perl-xml-simple=2.25=pl526_0
  - perl-yaml=1.27=pl526_0
  - prodigal=2.6.3=1
  - prokka=1.12=3
  - samtools=1.9=h8ee4bcc_1
  - spades=3.12.0=1
  - tbl2asn=25.3=0
  - boost=1.60.0=py36_0
  - bzip2=1.0.6=h14c3975_5
  - ca-certificates=2018.03.07=0
  - certifi=2018.10.15=py36_0
  - curl=7.62.0=hbc83047_0
  - expat=2.2.6=he6710b0_0
  - fontconfig=2.13.0=h9420a91_0
  - freetype=2.9.1=h8a8886c_1
  - giflib=5.1.4=h14c3975_1
  - icu=58.2=h9c2bf20_1
  - jpeg=9b=h024ee3a_2
  - libcurl=7.62.0=h20c2e04_0
  - libedit=3.1.20170329=h6b74fdf_2
  - libffi=3.2.1=hd88cf55_4
  - libgcc=7.2.0=h69d50b8_2
  - libgcc-ng=8.2.0=hdf63c60_1
  - libgd=2.2.5=hceca4fd_3
  - libpng=1.6.35=hbc83047_0
  - libssh2=1.8.0=h1ba5d50_4
  - libstdcxx-ng=8.2.0=hdf63c60_1
  - libtiff=4.0.9=he85c1e1_2
  - libuuid=1.0.3=h1bed415_2
  - libwebp=1.0.0=h222930b_1
  - libxml2=2.9.8=h26e45fe_1
  - ncurses=6.1=hf484d3e_0
  - openjdk=8.0.152=h46b5887_1
  - openssl=1.1.1=h7b6447c_0
  - perl=5.26.2=h14c3975_0
  - pip=18.1=py36_0
  - python=3.6.7=h0371630_0
  - readline=7.0=h7b6447c_5
  - setuptools=40.5.0=py36_0
  - sqlite=3.25.2=h7b6447c_0
  - tk=8.6.8=hbc83047_0
  - wheel=0.32.2=py36_0
  - xz=5.2.4=h14c3975_4
  - zlib=1.2.11=ha838bed_2

Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.
conda env create -f bio.yaml

I hope this gives you some help how to process your data. If you need some information just let me know ;)


Is it possible to merge scRNAseq data from experiments with different design?

It’s right there on the cellranger manual:

#aggregate results of counts for separate samples
cellranger aggr
#analyse the combined results
cellranger reanalyze

Note: I’m not sure this applies to the above case where batch effects could be an issue. However, this is how the analysis is performed for 10X data from multiple runs in principle.

cellranger aggr

cellranger reanalyze


Enrich for genes expressed in the plasma membrane

You can do functional enrichment analysis, for instance with gProfiler in R and see if your marker genes are enriched in any GO term related with the plasma membrane.

Alternatively, you can generate a list with genes whose product are located at the plasma membrane and perform Gene Set Enrichment Analysis, for instance with fgsea in R.


Adding external annotation to the phylogeny tree

there are nice features in R ggtree package. You may use the geom_cladelabel() function for your purposes. But I haven't find yet a way to automatize clade annotation. To do it, you have to be sure of the clade homogeneity before removing tiplabs. here are some commands to do clade label annotation:

annot<-read.delim("annota.txt", header = T)

ggtr<-as.treedata(pars.a)
g<-ggtree(ggtr)
#find the node number corresponding to the common ancestor of the top clade in you tree
MRCA(ggtr, "AAL40400", "AAF25519")
#here, check that all accessions belong to the same organism/group
subset(annot, label=="AAL40400" | label=="AAF25519")
#in this case, sequences are related to bovine.
g<-g + geom_cladelabel(node=89, label="Bovine", align=F, fontsize=4)
plot(g)

#another example
MRCA(ggtr, "AAA46456", "AAF97738")
subset(annot, label=="AAA46456" | label=="AAF97738")
g<-g+ geom_cladelabel(node=73, label="Rodent", align=F, fontsize=4)
plot(g)


Using preprocessing/alignment functions on the server

The short answer is to use conda. In the bioconda channel we have most of the tools used in bioinformatics, such as samtools. You do not need administrator permissions to install conda or packages with it, so your lack of sudo ability is not an issue.

  I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error

Installing packages on your local machine won't affect the remote server unless they share a file system, so this is expected.

As an aside, BAM files aren't really "raw", they're considered processed data (fastq/fast5 files are considered truly "raw" data in the sequencing world).


Question: How to simulate 100k samples having 40 million SNPs in a proportion of case:control=30:70?

First comment: you should be using real data, rather than simulated data, as we are not very good at modelling genetic variation and recombination at a genome-wide scale. A way that this can be done with GCTA is mentioned in the PLINK manual.

I'll do some calculations to estimate the size of data you need to store:

Assume that this plink-formatted file is in tped format (one row per SNP), and the entire genome is diploid. That means that a single SNP/individual pair will take up three bytes. You want 30,000 cases and 70,000 controls, so each line from each sample will take up 300kB (plus a bit of change for the header). With 1 thousand SNPs, that's 300MB per sample; with 1 million SNPs, that's 300GB per sample; with 10 million SNPs, that's 3TB per sample; with 40 million SNPs, that's 12 TB for one sample.

Reading the plink manual, I see that the BED format is a two-bit compression for each SNP/individual pair, so one byte per 4 individuals, or about 12x compression (i.e. 1TB for one sample).

Okay, that's a lot. It's well within a 256TB storage limit, but might hit file size limits depending on the file system.

If there is a data generation issue with PLINK, it might be worth trying (as you suggest) splitting the input for simulation then merging afterwards. My guess is that doing the simulation on fewer SNPs (rather than fewer individuals) would work better. Plink has a --bmerge option for merging the compressed binary file formats.

I usually approach problems like these by working my way up from something that I know works. Try first merging two datasets of 1000 SNPs. If that works, try larger datasets. Rinse and repeat until the computer crashes, you get bored, or you reach your goal.


How to get matching pattern along with ID in a single command in grep?

Use this Perl one-liner:
perl -lne '$id = $1 if /^>(.+)/; ($m) = /(CC[GT]AAA[GC][AC]TT[GC])/; print join "\t", $id, $m if $m;' input.fasta

The Perl one-liner uses these command line flags:
-e : Tells Perl to look for code in-line, instead of in a file.
-n : Loop over the input one line at a time, assigning it to $_ by default.
-l : Strip the input line separator ("\n" on *NIX by default) before executing the code in-line, and append it when printing.
SEE ALSO:
perldoc perlrun: how to execute the Perl interpreter: command line switches


Quantify low complexity regions in DNA sequences

"Complexity" is also known as entropy. A quick and dirty way is to use a file compression utility - zip will get good compression for a repetitive sequence.
You can also look at the k-mer distribution for some suitable k: in a repetitive sequence it will be biased. Here is a tool that give a local plot: https://rdrr.io/github/vsbuffalo/qrqc/man/kmerEntropyPlot.html


Presuming we consider "human readable", "easily parsable", and "quickly queryable" to be objectively good qualities (and if not, I worry for the future):

Text-based: It's absurdly common to want to use grep or awk on annotations. Sure, one could make variants of these that are binary-format aware, but why reinvent the wheel. Of course text files don't innately allow region-based queries of their contents, so on to point 2...

This should further explicitly be line-based. Fields would need to be tab-separated (ideally it'd use the ascii record separator character instead, but I fear that ship has long ago sailed).

A strictly defined line order: This is one of my personal pet-peeves about BED/GFF/GTF files. If you're going to make a text-based annotation format, everyone would end up ahead if said format made explicit that it should be sorted. This then allows things like tabix to be used so the "query a region" problem is then solved. But I would go further than that. The issue with things like GFF is that there are multiple inter-dependent lines and there's no strict rule about whether a parent line absolutely must come before a child line. This just makes implementing things a nightmare and more often than not a randomly ordered file will just break tools. Since GTF or GFF-style interdependent lines likely be how any annotation format works, this ordering between lines with the same start position should be made explicit in the format.


How can I run the fdnadist command from the EMBOSS package?

To find emboss packages type wossname, then hit enter twice.  You will see all the programs.  There is one called distmat.  Then do tfm distmat, to get the fkn manual for distmat.


pysam pileup: what reads appear in the pileup?

The two entries in the sam file represent mate pairs for the given ID.  You can tell the difference based on the sam flags.  Picard has a nice tool to determine the meaning of the flags.

You are missing coverage within the first read.  I think that it is related to the low quality bases.  samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.


Error in createMAE function: non-unique values when setting 'row.names' in TCGA LIHC Data

The error message is pretty clear: R requires unique row names on data frames and some of the "values" that you are trying to set as row names are non-unique.

You will need to provide a vector i) of length the same as your row number, ii) that is composed of non-unique values.


I would also recommend to take a look at pomegranate, a nice Python package for probabilistic graphical models. It includes solvers for HMMs and much more. Under the hood it uses Cythonised code, so it's also quite fast.


Question on nanopore sequencing data process pipeline (cDNA-PCR)

I recommend that you look at the most recent version of the protocol, which includes example adapter and barcode files:

https://www.protocols.io/view/stranded-mapping-from-long-reads-syheft6?step=2

[check in the protocol to see if there's a more recent version]

These files contains the public sequences that ONT have released for their sequencing kits (as briefly mentioned in the "Before Starting" section of the protocol). The sequences are available on the ONT community, but I have converted them into FASTA files for use in my protocol.

The adapter and barcode files were combined into a single FASTA file in the previous version of my protocol. For the PCR cDNA barcoding, it looked something like this:

>ONT_VNP
ACTTGCCTGTCGCTCTATCTTCTTTTT
>ONT_SSP
TTTCTGTTGGTGCTGATATTGCTGCCATTACGGCCGGG
>BC01
AAGAAAGTTGTCGGTGTCTTTGTG
>BC02
TCGATTCCGTTTGTAGTCGTCTGT
>BC03
GAGTCTTGTGTCCCAGTTACCAGG
>BC04
TTCGGATTCTATCGTGTTTCCCTA
>BC05
CTTGTCCAGGGTTTGTGTAACCTT
>BC06
TTCTCGCAAAGGCAGAAAGTAGTC
>BC07
GTGTTACCGTGGGAATGAATCCTT
>BC08
TTCAGGGAACAAACCAAGTTACGT
>BC09
AACTAGGCACAGCGAGTCTTGGTT
>BC10
AAGCGTTGAAACCTTTGTCCTCTC
>BC11
GTTTCATCTATCGGAGGGAATGGA
>BC12
CAGGTAGAAAGAAGCAGAATCGGA
>RB12A
GTTGAGTTACAAAGCACCGATCAG

I separate out adapters and barcodes in the newer versions of the protocol because I was getting too many incorrect hits.


How to read and interpret a gene expression quantification file?

The first column contains Ensembl gene identifiers, and the suffix is a version number that can be used to track changes to the gene annotations over time. From the Ensembl Stable IDs documentation:

  Ensembl annotation uses a system of stable IDs that have prefixes based on the species scientific name plus the feature type, followed by a series of digits and a version e.g. ENSG00000139618.1. The version may be omitted.

Following the first link you provided leads to a page with details for the file 2edcaaa7-63b4-40b4-abbe-5d7a84012e60.FPKM-UQ.txt.gz. The first thing that caught my eye about this filename was FPKM, or "fragments per kilobase of exon per million reads", which is a commonly used unit of RNA expression. Since these are not raw read counts, there's no expectation that these values should be integers.

The best explanation I've seen of FPKM comes from a blog post written by Harold Pimentel of kallisto and sleuth fame. From the blog post:

  The interpretation of FPKM is as follows: if you were to sequence this pool of RNA again, you expect to see FPKM_i fragments for each thousand bases in the feature for every N/10^6 fragments you’ve sequenced. It’s basically just the rate of fragments per base multiplied by a big number (proportional to the number of fragments you sequenced) to make it more convenient.

More generally though, even when FPKM is not the unit of expression abundance used, most quantification methods and associated units of expression will not produce integer estimates.


What is the significance of the effect size on WOMAC scores for TKA?

Hedge's g measures the difference between two groups here experimental and  control in terms of combined (weighted) standard deviation. A higher g implies greater difference and lower g means lesser difference between  two groups.
It depends  on the study how the scoring of questions is proceeded with : the effect-size will be same. For example, answer to question item may be scored 1 to 10 depending on whether question item is written in positive  sense (the apple is good for health) or negative sense (the apple is not good for health). Positive Hedge's g does not indicate larger effect for the treatment group. THE negative Hedges' g for the other two studies indicates positive effect for the treatment group if the difference has been computed by  deducting experimental group mean from the control group mean (which is less than exp. Group Mean). Overall, I find that a knowledge of questionnaire scaling and scoring of scales and question items could be very helpful in general and in doing a meta-analysis effectively.


First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway.

Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the phred scaled probability of the alignment/base being wrong. It often amuses me that we question mapQ but take baseQ for granted. BaseQ is also scaled and discretized differently; even fewer people know how Illumina/pacbio/nanopore/historical sequencers estimates baseQ.

On the second question, Section 2 of the MAQ supplementary explains the theoretical aspects of mapQ, which is still correct today. Briefly, mapping quality consists of three components: 1) the probability of contamination; 2) the effect of mapping heuristics and 3) the error due to the repetitiveness of the reference. Only 3) can be modeled theoretically.

In case of bwa-mem, if we assume the matching score is 1, type-3 error is estimated with:

$$10/\log10\cdot[\log4\cdot(S_1-S_2)-\log n_{\rm sub}]$$

where $S_1$ is the best alignment score, $S_2$ is the second best and $n_{\rm sub}$ is the number of suboptimal alignments. Factor $\log 4$ comes from the scoring matrix. Factor $10/\log10$ is the Phred scale. This equation assumes gap-free alignment and is very close to Section 2.5.2 in the MAQ supplementary. It is ok-ish for short reads, but often overestimates for long reads. I am not aware of a practical approach in general cases. In addition to this method, you can estimate mapQ by read simulation: just try to find a function such that it fits empirical mapQ. Some have tried machine learning, too.


Convert CD (cluster of differentiation) into gene symbol

I ended up using option 4 and parsing the columns for "Synonym" and "Previous Symbol".

https://www.genenames.org/cgi-bin/genefamilies/set/471


Master of Pores Installation Mac OS Mojave 10.14

you should not move nextflow within ~/master_of_pores/NanoPreProcess/bin .
Just add nextflow in your path.
export PATH=$PATH:where_is_your_nf_executable

Feel free to open an issue on github if you need more help.


Purpose of correlating docking score with catalytic efficiency(Kcat/km)

What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?

This is specific to the enzymology of paper you are citing. If the solubility of your substrate is lower than your Michaelis constant (K_M) you cannot fit to a Michaelis-Menten curve and can only get the slope of the initial straight part of the curve, which is the catalytic efficiency (k_cat/K_M). Ideally K_M would be better. But in generally catalytic efficiency ought to correlate with well calculated ddGs of ligand binding.

  I saw some papers are doing this, is it to validate the docking protocol?

Totally. Docking finds were ligands bind. Although the energetics don't match too well the experimental kinetics as of 2019.
The catalytic efficiency (k_cat/K_M) is the ratio of the turnover (k_cat) over Michaelis constant (K_M). The Michealis constant is in turn the ratio of on-rate over off-rate plus turnover. The ΔΔGibbs of the ligand ought to correlate with the ratio of the on-rate/off-rate, while the turnover is unrelated.
The turnover should correlate with the difference in energy between the substrate bound form and the transition state. This is the realm of theozyme design and requires correctly protonated ligands and protein, but when it is you get a nice transition state hill (e.g. personal example).

There are a few "ought to" in the above. Docking and force-field energy scores are not precise. The paper you refer to is doing very advance docking (akin to dynamic undocking) and, not accidentally, the ligands are very flexible and rather soluble.

  Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

A positive correlation coefficient just means that generally both the catalytic efficiency and docking scores increases or decrease together. If a metric "worked backwards", such as ΔΔGibbs (more negative is better), you would want a negative correlation. What you want is a strong correlation, with the sign a result of the nature of your metrics.


Picard, like most tools, prints logs to stderr.
It should be in the Galaxy history. Are you using a version of picard from the toolshed (if so, which one)? You might set cleanup_job to never in the config, which will prevent galaxy from cleaning the working directory. This will enable seeing what the job actually produced.

Most likely there's more in the log file that you're not showing is us and picard itself never ran to completion. In lieu of that information, what I wrote above is usually sufficient to debug things.


Adding legend to Botton and top annotation in complexheatmap

Answer from @ram-rs, converted from comment:
For the Group, try using it without the anno_simple. Just Top_ha = HeatmapAnnotation(Group=top$Group, "Total Mutation Burden" = ...), nothing else for Group. See if that generates a legend at all. We can customize later.hing else for Group. See if that generates a legend at all. We can customize later.


How to validate that BAMs have been downloaded correctly?

samtools quickcheck is all you need. From the manual:

Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only).
Data in the middle of the file is not read since that would be much more time consuming, so please note that this command will not detect internal corruption, but is useful for testing that files are not truncated before performing more intensive tasks on them.
This command will exit with a non-zero exit code if any input files don't have a valid header or are missing an EOF block. Otherwise it will exit successfully (with a zero exit code).


Bedtools get fasta and ORF from a blastX run

Blast uses a 1-based coordinate system. Bed is 0-based half open. Your code needs to convert appropriately. See for example: https://www.biostars.org/p/317112/ .


I have previously estimated tumour purity with the EXPANDS an inferred tumour heterogeneity program which is designed to calculate the number of clonal subpopulations in matched tumour/normal samples.  The purity is essentially the size of the largest subpopulation identified in that sample - this is discussed in the programs FAQ.  In addition to a matched tumour / normal exome or genome sequencing you'll also need match somatic copy number as a seg file as input.  Other programs also exist for this sort of inferred heterogeneity analysis - some of which may also give you a measure of purity are: Absolute, ThetA, SciClone, CHAT, PyClone and Canopy. A more complete list looks to be here.

The only other thing I'd suggest is have an estimate via some orthogonal measure of the purity to try at least judge how the different methods are performing with you data, e.g. you might have clonality indications from cytogenetics/FISH or histology work, or perhaps FACS.  Picking samples known to be pure / impure via one or more of these might help you get a handle on how well various methods are performing.

As regards estimation without a normal sample the only program which looks like it might help is QuantumClone, but that requires more than one tumour sample form said patient either spatially or temporally distinct to perform the analysis.

Also If you have low coverage DNA sequencing CNAnorm looks like it could be useful with just a single sample.


The positions you see in the x-axis of this plot are positions on your reads, not on the reference genome. Since each read comes from a random position on the genome, the frequency of A/T/C/G on each base position should reflect the base composition of the entire genome. Therefore, the percentage of each base should be approximately constant across the read.

Hope that helps!


Could it be that these are reads with their mates found on other chromosomes? That's the part of the IGV documentation which seems to make sense in this context:

for paired end reads that are coded by the chromosome on which their mates can be found


Combine Fastq by writeFastq is not working properly

You can drastically simplify your effort because FASTQ files are simple text files; therefore, standard UNIX test file tools work:

cat refs/*.fastq > combined.fastq

And this even works with .fastq.gz files — no unzipping required.

Apart from this, your R code contains several errors. For instance, your variable names change throughout, and you’re using the fsep parameter wrongly. In fact it’s unclear to me what the code is supposed to do. At an approximation (and untested), the following should work:

fastq_dir = file.path(getwd(), "refs")
fastq_files = dir(fastq_dir, pattern = "\\.fastq$", full.names = TRUE)
outfile = file.path(dirname(fastq_dir), "combined.fastq")

for (file in fastq_files) {
    fq = readFastq(file)
    writeFastq(fq, outfile, mode = "a")
}

This is equivalent to running the shell command shown above, just in a more convoluted way (parsing and re-generating the FASTQ files).


What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve?

Typically I use the Mauve viewer to visualize the alignments.  Are there things you find lacking in that viewer?

There are multiple threads about converting xmfa to other formats


Error in quantile.default(x, p = p) in EgdeR calcNormFactors

Ok, it turned out that some rows and/or columns contained all values as zeros and that was the reason for the error. I wrote a function to fix that:

def fix_rawcounts(rawcounts, cell_ids, gene_symbols):
    '''The function checks whether there are
    any rows or columns with all zeros
    and removes them'''
    print('Checking whether rows or columns contain all zeros')

    # Checking columns
    count_arr = np.apply_along_axis(np.count_nonzero, axis=0, arr=rawcounts )
    rawcounts = rawcounts[:,count_arr != 0]

    # Checking rows
    count_arr = np.apply_along_axis(np.count_nonzero, axis=1, arr=rawcounts )
    rawcounts = rawcounts[count_arr != 0,:]

    print('finished checking zero rows and columns')
    return rawcounts


Is there a tool for synonymous recoding of DNA sequences?

There seems to be a bioconductor package named synmut which does just that.
Apparently, it can take codon usage into account when generating synonymous mutations.

Copy pasting an example that matches your use case from the official doc:

library("SynMut")
filepath.fasta <- system.file("extdata", "example.fasta", package = "SynMut")
filepath.csv <- system.file("extdata", "target_regions.csv", package = "SynMut")
region <- read.csv(filepath.csv)
rgd.seq <- input_seq(filepath.fasta, region)
mut.seq <- codon_random(rgd.seq, keep = TRUE)
Biostrings::writeXStringSet(get_dna(mut.seq), "rgd.fasta")


Which gene I should select from this qqplot

Given that TP53 is the most significant and is already known to have driver mutations in cancer it would seem to be the logical choice.


How to identify unknow bacterium species from whole genome genetic sequence

You can use KmerFinder and Kraken2 for taxonomy assignment


Is removing samples based on clustering for downstream analysis a right choice?

I would be very hesitant to blindly exclude those samples based on the clustering. Check to see if the clusters actually denote some sort of batch effect, since it's not like all of the TCGA datasets were processed at the same place or time. Check for clinical covariates too. Also, do a PCA and check the genes with high loading on the relevant PC. If it appears that the genes are indicators for cancerous state, then it's likely the samples are just mislabeled and you can reasonably exclude them. Otherwise keep them in so your results better represent the underlying biological variability.

BTW, your heatmap and MDS plots are what I'd expect from a cancer dataset. You occasionally see a tumour with the normal samples and vice versa, since the problematic genes aren't always in the top 10% and there is likely some heterogeneity in the cancer samples.


Identifying nested/overlapping reading frames in nucleotide sequences

An RNA sequence has six possible reading frames depending on the offset (0, 1, 2) and the direction (forward, reverse complement). All frames can potentially be used for  protein translation, but only one frame per protein (i.e. there is no shifting going on when translation is happening, excluding shifts introduced by splicing).
If you need to find all the possible protein sequences within an RNA sequence, there are plenty of programs already available (e.g. getorf and transseq from the EMBOSS tool kit).


Tools for creating gene coexpression networks

After looking around, it seems like the Camoco library is suitable.


genes with highest RPKM/FPKM for a human RNA-seq experiment?

As asked as it is, the answer is probably no.
Indeed, the highest expressed RPKM/FPKMs will be different from one condition (or tissue) to another one for example.
Then, you may also have technical artefacts due to the wet-lab part or to the normalisation. For example, mitochondrial genes are often reported in the top expressed genes.
Now, if you want to compare your results to an outside source, I recommend checking the EBI gene expression Atlas* and more specifically the Baseline expression experiments.
Pro-tips:

You can filter the genes based on their nature (eg, mRNA) and/or their expression.
You can also visualise the overall expression of a specific gene across several experiments.

*Disclaimer: I am related to this project myself.


sva for RNA-Seq data without known phenotype

You can just remove the first component from the dataset by setting the first eigenvalue to 0 in the diagonal matrix and then multiplying the SVD matrices.

I am not sure which R code you are using to estimate the components (I have a comparison of PCA code in R on https://github.com/aedin/ODSC_2018/tree/master/PCA_vignette)

In ade4 R package, there is a function reconst() which reconstructs data from a selection of components: https://rdrr.io/rforge/ade4/man/reconst.html

Happy to provide code to do this in other PCA versions if needed
Aedin Culhane


How to indicate the END of a haplotype block in VCF?

According to the VCF specification, haplotype blocks are defined by the "Phase Set" [PS] tag:

  PS : phase set. A phase set is defined as a set of phased genotypes to which this genotype belongs. Phased genotypes for an individual that are on the same chromosome and have the same PS value are in the same phased set. A phase set specifies multi-marker haplotypes for the phased genotypes in the set. All phased genotypes that do not contain a PS subfield are assumed to belong to the same phased set. If the genotype in the GT field is unphased, the corresponding PS field is ignored. The recommended convention is to use the position of the first variant in the set as the PS identifier (although this is not required). (Non-negative 32-bit Integer)

A good example of this can be found on a 10X Genomics support page about phasing:

CHROM   Pos REF ALT GT  PS
chr1    1000    A   C   0|1 1000
chr1    1010    T   G   1|0 1000
chr1    2000    C   T   0|1 2000
chr1    2005    T   G   0/1 2000
chr1    2008    G   C   0|1 2000

  In this example we have two phase blocks, denoted by PS=1000 and PS=2000. PS=1000 spans position 1000-1010, and PS=2000 spans position 2000-2008. In PS=1000, haplotype 1 contains the REF A allele at position 1000, and the ALT G allele as position 1010, while haplotype 2 contains the ALT C allele at position 1000 and the REF T allele at position 1010.

Note that this phase set identifier needs to be a single numeric value, which means that sub-blocks and overlapping blocks cannot be defined by this field.


What does instantaneous rate matrix mean?

Reversible At its simplist the GTR rate, is the General Time Reversible model and most importantly infers the matrix is "reversible", thus as many A mutations will go to C as C mutations go to A and the number of mutations is relative to the G-T mutation rate. The G-T mutation rate is one of the least common mutations (transversion) so all other values are above 1. If you plotted this as a matrix the top right-hand matrix is identical to the bottom left hand matrix because its reversible, i.e. A <-> C . This is cool because you don't need to tell the algorithm which sequence is the root.
Directional The "instantaneous" matrix is a directional matrix, thus A -> C is not the same as C -> A. In order to obtain an accurate directional matrix you need to assign a root, its like lets start from the deepest node in the tree (the one between the root/outgroup and all ingroup sequences) and now - if that is correct - we can remove the restrictions of reversibility and impose directionality on the mutation process. This is needed if you are performing molecular clock type calculations. However, the key is assigning the root/outgroup - if it assigns it e.g. by mid-point rooting that is not necessarily a good idea.
Advice. My advice is if you are building a basic tree and haven't assigned a root/outgroup, obviously using GTR which is a good move, then you don't to worry about the intananeous matrix because it might not be meaningful. If you are running a molecular clock/dating/rate estimate type model and a very clear what the root/outgroup is and that is specified in the input, then the instananeous matrix becomes important.

The mutation rates of GTR and instananeous are not directly comparable because the GTR rate is with respect to G <-> T mutation, whilst instananeous is primarily concerned with the directionality of mutation and is a 12-parameter model, as opposed to GTR which 6-perameters, where one of those parameters is the reference.
Thus G-T is reversible in GTR (the reference) and so e.g. A<->C is reversible with respect to G<->T.
In instananeous matrix G -> T and T -> G are two different rates as are A-> C and C -> A, so one approximation based on an assumption with respect to a second approximation based on the same assumption can't directly equate to a mutation rate which removes that assumption.


Truncating branch length values of Phylogenetic tree with biopython

Apologies, try the following,

allclades = list(pars_tree.find_clades(order='level'))
for i, blength in enumerate(allclades):
        printf(blength=%.2f, blength)

"allclades" is a dictionary rather than a simple array.

If this fails I would examine the allclades data structure either using pprint (directly) or using printStruct via the following code,

>>> d = [{'a1':1, 'a2':2, 'a3':3}, [1,2,3], [{'b1':1, 'b2':2}, {'c1':1}], 'd1', 'd2', 'd3']
>>> printStruct(d)

def printStruct(struc, indent=0):
  if isinstance(struc, dict):
    print '  '*indent+'{'
    for key,val in struc.iteritems():
      if isinstance(val, (dict, list, tuple)):
        print '  '*(indent+1) + str(key) + '=> '
        printStruct(val, indent+2)
      else:
        print '  '*(indent+1) + str(key) + '=> ' + str(val)
    print '  '*indent+'}'
  elif isinstance(struc, list):
    print '  '*indent + '['
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ']'
  elif isinstance(struc, tuple):
    print '  '*indent + '('
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ')'
  else: print '  '*indent + str(struc)

Once you know the structure you unravel it accordingly to get to the branch lengths and then issue a printf command to round them to 2 dp.

Goodluck


You can exclude the which statement and get the entirety of all of the chromosomes. Note however that R is a very bad platform to use for this, since you end up reading everything into memory (I hope the BAM file isn't huge).
Note that you're seriously reinventing the wheel. There are a number of packages that will directly produce the plot you want, such as bamPEFragmentSize from deepTools.


The cause of the error is that your input file is truncated or /tmp is running out of space. If you can do samtools view -H input.bam without error (reading the header also checks for the magic number at the end of the file) then it's like the that /tmp is filling up.


swap genotypes for specific fields of a vcf file

If we can always assume your fields will be GT:AD:DP:GQ:PL, you can use this relatively simple, brute approach awk script:

awk -vOFS="\t" -F"\t" '{
  ## just print the header lines
  if(/^#/){
    print;
  }
  else{
    ## The number of samples will be the extra number of fields
    ## in addition to the default VCF 9
    sampleNum = NF - 9
    ## Fix each sample
    for(i=1; i<=sampleNum; i++){
      ## get the fields
      num=split($(9+i), fields,":");
      ## genotype
      if(fields[1] ~ "0/0"){ fields[1] = "1/1"}
      else if(fields[1] == "1/1"){ fields[1] = "0/0" }
      else if(fields[1] == "0/1"){ fields[1] = "1/0" }
      else if(fields[1] == "1/0"){ fields[1] = "0/1" }
      ## AD
      split(fields[2],values,",")
      fields[2] = values[2]","values[1]
      ## PL
      split(fields[2],values,",")
      ## Now set these as the sample data
      $(9+i) = fields[1];
      for(k=2;k<=num;k++){
        $(9+i) = $(9+i)":"fields[k]
      }
    }
    print;
  }
}' file.vcf

On your example input, this returns:

Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:0,48,720 1/0:8,14:22:99:167,0,347    .....


Plotting Ramachandran Plot from more than one PDB file

You could probably use PyRAMA for that. The only problem is that by PyRAMA plots to display, while you probably would like to save to a file.

One way to get around this is to:

Patch PyRAMA to allow ramachandran function to accept output image filename.

def ramachandran(file_name_list, output_image=None):
    ...
    plt.tight_layout()
    if output_image is not None:
        plt.savefig(output_file, dpi=300)
    else:
        plt.show()

If this works as expected I would really encourage you to submit a pull-request (or at least open an issue) on github.
Now you can write your main script:

from glob import iglob

from rama_lib.ramachandran_calc import ramachandran

for pdb_filename in iglob('/my/dir/with/pdb/', '*.pdb'):
    img_filename = pdb_filename.rsplit('.', 1)[0] + '.png'
    ramachandran((pdb_filename,), img_filename)

I did not test the code, but I hope you get an idea.

There are also other ways to make a ramachandran plot in python. Please check: 1, 2, 3.


Is there a Python package to convert InChi to molecular structures?

Yes, RDKit can be used, however, if you installed it with conda it will not work out of the box for inChi key fetching. You can either spend some time installing the missing bit or use something else. It is not like Mol2 support in RDKit, which shouldn't be used, it is just not worth the bother.

I would suggest PubChem's RESTFul API.

Say

 import requests
 ethanol = 'LFQSCWFLJHTTHZ-UHFFFAOYSA-N'
 r = requests.get(f'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{ethanol}/property/CanonicalSMILES/JSON').json()
 smiles = r['PropertyTable']['Properties'][0]['CanonicalSMILES']


Which is the amino acids classification when analyzing the solvent accessibility of a protein?

There is some confusion in the question.

Solvent accessibility of a residue, i.e. the area accessible to the solvent (water), depends on the 3D structure of the protein. It's not a property of  an amino-acid (chemical compound such as alanine, leucine, ...).

There may be a correlation -- amino-acids with hydrophobic side chains are more likely to be buried, but in general amino-acids (chemical compounds) cannot be classified as buried or exposed.

I don't have access to the first article, but I tried running the prediction server described in the second linked article, and I got a result

GSGMQIFVKTFRDRLRNLLPQTITLEVEPSDTIENVKAKIQDKEGIPPDQQVLIFSRKRL
CCCCEEEEEECEEEEEECCCCEEEEEECCCCCHHHHHHHHHHHHCCCHHHEEEEECCEEC
EEEBeBbBebbEbEbEebEEeEbeBEBeeEeeBeeBbEeBeEeEEbEeEeBeBbbEEeeb

...

Porter (Secondary Structure):
H = Helix   (DSSP classes H, G and I)
E = Strand  (DSSP classes E and B)
C = Coil    (DSSP classes S, T and .)

PaleAle (Relative Solvent Accessibility):
B = very buried      (<=4% accessible)
b = somewhat buried  (>4% and <=25% accessible)
e = somewhat exposed (>25% and <=50% accessible)
E = very exposed     (>50% accessible)

In one place in the sequence you can see QQ, predicted as eB.
So clearly glutamine is NOT assigned to a particular RSA class.

Dividing RSA (relative solvent accessibility, 0-100%) into categories is arbitrary, I don't think there is any widely used convention. The easiest way would be to just use numeric values.


How to generate the json files for Cromwell workflow execution?

Unfortunately generating inputs for Cromwell is not something that is handled by Cromwell itself or by GATK, and generally needs to be scripted to fit your custom needs.

Despite this being a bit arduous, there is a very good reason for this. WDL itself is a highly polymorphic workflow language, where inputs can be vastly different between each workflow. Cromwell ONLY understands how to parse the workflow and the inputs, but does not know the context of what it is running (ie this is a paired mapping job), nor does GATK know that it is running within a WDL.

Additionally, I've encountered many different cases where there are not just 2 fastq files but, due to older sequencing technologies, 10 or more for the same sample. Unless GATK/Cromwell knows how your specific samples are arranged there is no easy way for them to generate your inputs.

There are tools that would possibly help you simplify scripting better. WOMTool developed by the Broad has an inputs command which will take an  arbitrary workflow and generate a JSON template of the inputs, which you can then use in your own scripts.

As an example, imagine you have this WDL workflow in a file called my.wdl:

version 1.0
workflow foo {
  input {
    Array[Pair[File,File]] fastqFiles
  }
}

When you run WOMtool's inputs command, it produces a template JSON listing all inputs and specifying the variable types expected for each:

$ java -jar womtool.jar inputs my.wdl
{
  "foo.fastqFiles": "Array[Pair[File, File]]"
}

You can use this template to produce actual JSONs for each of your samples.

Note that you can also get the template written to a file as follows:

$ java -jar womtool.jar inputs my.wdl > my.inputs.json


EGA remote download via snakemake fails due to EGA_CLIENT_ID specification

Okay, I finally solved it. The reason was the EGA_CLIENT_ID variable. In the pyega3.py script, there is a definition for the client_id for pyega3 v3.0.44 ("f20cd2d3-682a-4568-a53e-4262ef54c8f4"). After setting the EGA_CLIENT_ID to this value login was successful.

os.environ["EGA_CLIENT_ID"] = "f20cd2d3-682a-4568-a53e-4262ef54c8f4"

I think this is a bug and I filed it on Snakemake's git repository - the official documentation just requires definition of the EGA_USERNAME and EGA_PASSWORD variables, and while definition of EGA_CLIENT_SECRET will probably need to be done by the user, the client ID should be pulled automatically from the pyega3 version installed locally.

Of note, after successful login, the example download from the documentation failed with the following:

Missing input files for rule a:
ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai

However, downloading files that I have access to worked, indicating that the paths to the files from the documentation are wrong or are access-protected.


How to output only the overlapping region of two files?

bedtools is my go-to program for operations on genomic intervals. In particular, the bedtools intersect operation is what you're looking for here.

$ cat fileA
chr1    25      50
chr1    75      200
$ cat fileB
chr1    10      60
chr1    80      90
$ bedtools intersect -a fileA -b fileB
chr1    25      50
chr1    80      90

Make sure that the values are separated by tabs rather than spaces, or bedtools will have trouble reading the files.


Is there a way to do GWAS on phenotype data that is not normally distributed?

Yes, that can be done.
A common approach for these types of analysis is to carry out a transformation on the data first in order for it to have a normal distribution (the most general approach called an inverse normal transformation, among other names), then run the test on the transformed data.
As an example, for data that are positive and zero-skewed, a log transformation may be appropriate.
The transformation approach does require variation in the data values, though. For a trait that is exactly zero for most values, it might be better to treat it as a binary/categorical trait for the purpose of association testing, rather than a quantitative trait. Most association testing packages should allow associations for binary traits to be tested.


It looks like the image you downloaded is missing the partition with the label rostlab-data. Your blkid output shows that you have the local-storage partition, but not the rostlab-data which is automatically added to your /etc/fstab.

I am afraid you'll have to contact the developers and let them know. Sounds like a bug.


You are looking at the supplementary data of a paper. That seems to have given you a list of features, and some information about those features. Specifically, you seem to have a list of two types of element:

"peg". Based on the information there, I assume "peg" stands for "protein encoding gene". Note that all peg lines have a peptide sequence (not seen in your question but present in the original link) and their "function" column is always some sort of protein function (e.g. "hypothetical protein" or "Transposase").

"repeat". These are presumably repetitive elements.

For each line, you are given an identifier for the sequence in which this element can be found. For example, your first line is on gi|1045318032|gb|MAZB01000333.1|. That is a GI (gi|1045318032) and a GenBank identifier (gb|MAZB01000333.1). If you go to the NCBI website, and search for 1045318032, you will find that it is the Candidatus Erwinia dacicola isolate Erw_SC contig_333, whole genome shotgun sequence.
So, your first line is telling you that the authors have identified a protein encoding gene on sequence gi|1045318032 and the gene starts at position 511 of that sequence and ends at position 2. This might sound odd, but that's because the gene is on the negative (-) strand and positions are calculated with respect to the positive strand. So, if you were to count on the negative, the gene would start at 2 and end at 511.
The same goes for all the other lines you have there. Just look up the GI number and you can know what sequences they are referring to.
Finally, the second column, the "feature_id" looks like something specific to this figshare service that is hosting the data. I dont' know what it is, but I would guess it's an identifier for this particular feature. I don't think it has any relevance outside http://figshare.com.


Identifying proteins in a concatenated protein sequences

The only way I can think is to break the protein apart on methionine residues, blast each amino acid locus and concatenate any locus that does not score a query protein position of '1'. The theory is simple every protein starts with a methionine residue, however there are plenty of M residues within the protein that don't represent the beginning of the gene.

Breaking the protein apart on M would be easy, but then automating and parsing the blast would be fiddly. Obvioulsy you will need to code it and if its a one off would it be worth it?


Obtaining all protein sequences with a particular domain architecture from Pfam

An alternative way is using the Hmmer website, and do a hmmsearch with the Pfam accession, forex. (PF02910). Navigate to the domain tab and find the domain architecture that you are interested in and press view scores. That way, only the sequences you are interested in will be downloadable under the tab downloads.

Step-by-Step guide:

1- Go to http://www.ebi.ac.uk/Tools/hmmer/

2- Search with hmmsearch http://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch

3- Use Accession search and type in PF02910

4- Navigate to Domain

5- To the right side of
9742 sequences with domain architecture: FAD_binding_2, Succ_DH_flav_C

Press View scores.

6- It should show (Your results have been filtered) and Query matches (9742)

7- Navigate to Download and download the sequences in the format of interest.

8- Repeat for domain architectures of interest.

I hope you find this helpful.

Kind Regards,

Sara


Encoding a 15bp DNA sequence into a shorter than 15 digits number

What you are asking is called a "hash" (wiki).
In python, get your oligo as a string and use the standard library method hash(oligo), this is a number that gets used for a lot of operations. It is however always a 18+1 digit number in decimal, so a different solution is needed.
There are four nucleobases, so binary conversion from scratch is easy and best option and each solution will be unique and small. Each base takes up 2 bits and you concatenate these and convert them to hex or decimal —the length should be 10 decimal digits. One issue is when you have a oligo of a different length which start with adenosine or not.
One issue with this is the order. I am not sure where it comes from but the bases are nearly often expressed in the order ATGC. So...
mapping = {'A': '00', 'T': '01', 'G': '10', 'C': '11'}

It is satisfying how the first digit is the W/S, while the second is R/Y.
Now simply covert and add them compressed = functools.reduce(lambda acc,base : acc+mapping[base], oligo, '') and covert the string to decimal like so int(compressed, 2).


I would guess that the real unit for measuring protein expression would be pM (pico-moles) or a count of protein molecules.

I'm not sure we have any methods to calculate such quantities though.


Which module to select for Pathway analysis based on Module trait correlation and pvalue?

If you can run pathway analysis programmatically on all modules, yes, do it on all modules. Otherwise, select the modules with strongest associations with each subtype, both positive and negative. I would keep the modules separate because they may have very different biological functions.


NCBI blast for exact match of a short sequence

UPDATE: while procrastinating, I worked up the below throwaway code into a more generally applicable (but still bare-bones) script.

ORIGINAL:
I feel like you might want to use something other than BLAST to do this.
BLAST is not really designed to solve this problem, I think it's just easy to use. Also, your N-containing sequence (16nt) is shorter than the version without Ns (21 nt), which can't be correct if you care about "exact" matches.
As a side note, the idea of exact matches is a bit tricky with Ns involved. The whole point of ambiguity codes is that things won't be exact.
And all of this is on top of the short-string problem- most alignment programs won't even return results that don't meet some minimum alignment score, which is exacerbated a bit by the Ns. For how BLAST handles ambiguity codes, you can see section 7.2 here.
Why not use regular expressions instead to search your sequence against the transcriptome?
# python
import re
from Bio import SeqIO

fwd = r'ATTG[ACGT]{4}GCAAACCA'
rev = r'TGGTTTGC[ACGT]{4}CAAT'  # reverse complement

# reference file
seqfile = SeqIO.parse("your_reference_file.fa", "fasta")  # a transcriptome reference fasta

for seq in seqfile:
    # maybe use re.match() to get coords but easier to show it with findall
    a = re.findall(fwd, str(seq.seq))
    b = re.findall(rev, str(seq.seq))
    if len(a) > 0:
        print(seq.id, a)
    if len(b) > 0:
        print(seq.id, b)

When I put this in a file and run it against a test case with ATTGAAAAGCAAACCA, it works:
$ python exact_match.py
the_contig ['ATTGAAAAGCAAACCA']

Alternately, you could just blast the 16 non-N versions of that sequence and combine them, as you've shown.


Figured it out! For anyone who finds this thread in the future and is wondering what is going on: HOMER, when using the -size given option, normalizes the y-axis to the number of basepairs in the intervals -- in this case, the length of the gene body. To get the "raw" values, as I wanted, you need to multiply by the length of the gene body to "undo" the normalization.


Reconstruct a protein's 3D structure from a pair-wise distance matrix

Finding coordinates of N points that preserve (as well as possible) pair-wise distances from an NxN matrix is called multidimensional scaling.

But I don't enough about the context of your problem to tell if/how it is applicable.


How to select a cutoff for interaction confidence in STRINGdb?

This depends on what you are trying to do and whether you value specificity over sensitivity. We can't tell you since it is entirely dependent on the biological question you want to answer.
However, I would recommend two things:

Don't use STRING. The creators of STRING made the choice to value sensitivity over all else, so they include any interaction they can get their hands on. That is sometimes fine, depending on what you want to do, but is more often a problem. Personally, I tend to avoid STRING as much as possible and instead use more curated databases like APID or IntAct. You can find a list of many PPI databases in the EBI's Psiqcuick View page. You might also be interested in my answer here which gives an example script using APID which can easily be modified to query other DBs.

Whatever you use, I recommend you filter by detection method. There are many different interaction detection methods ranging from those that only identify direct, binary interactions between two proteins (e.g. Yeast two-hybrid), those that also find proteins that are in the same complex even if there is no direct interaction between them (e.g. ChIP), to various non-experimental methods which are used to infer interactions.
Luckily, all of these have been collected into a controlled vocabulary which you can see in the Interaction Detection Method Ontology Lookup Service page. Go there, collect the interaction detection methods you consider good enough for whatever you want to do and filter your interactions using those.

For example, for a high quality human interactome I built for a project I used to work on, I was only interested in direct binary interactions, so I only kept interactions detected by the following methods:
    "MI:0008" => "array technology",
    "MI:0009" => "bacterial display",
    "MI:0010" => "beta galactosidase complementation",
    "MI:0011" => "beta lactamase complementation",
    "MI:0012" => "bioluminescence resonance energy transfer",
    "MI:0013" => "biophysical",
    "MI:0014" => "adenylate cyclase complementation",
    "MI:0016" => "circular dichroism",
    "MI:0017" => "classical fluorescence spectroscopy",
    "MI:0018" => "two hybrid",
    "MI:0020" => "transmission electron microscopy",
    "MI:0030" => "cross-linking study",
    "MI:0031" => "protein cross-linking with a bifunctional reagent",
    "MI:0034" => "display technology",
    "MI:0040" => "electron microscopy",
    "MI:0041" => "electron nuclear double resonance",
    "MI:0042" => "electron paramagnetic resonance",
    "MI:0043" => "electron resonance",
    "MI:0047" => "far western blotting",
    "MI:0048" => "filamentous phage display",
    "MI:0049" => "filter binding",
    "MI:0051" => "fluorescence technology",
    "MI:0052" => "fluorescence correlation spectroscopy",
    "MI:0053" => "fluorescence polarization spectroscopy",
    "MI:0055" => "fluorescent resonance energy transfer",
    "MI:0065" => "isothermal titration calorimetry",
    "MI:0066" => "lambda phage display",
    "MI:0073" => "mrna display",
    "MI:0081" => "peptide array",
    "MI:0084" => "phage display",
    "MI:0089" => "protein array",
    "MI:0090" => "protein complementation assay",
    "MI:0091" => "chromatography technology",
    "MI:0092" => "protein in situ array",
    "MI:0095" => "proteinchip(r) on a surface-enhanced laser desorption/ionization",
    "MI:0097" => "reverse ras recruitment system",
    "MI:0098" => "ribosome display",
    "MI:0099" => "scintillation proximity assay",
    "MI:0107" => "surface plasmon resonance",
    "MI:0108" => "t7 phage display",
    "MI:0111" => "dihydrofolate reductase reconstruction",
    "MI:0112" => "ubiquitin reconstruction",
    "MI:0114" => "x-ray crystallography",
    "MI:0115" => "yeast display",
    "MI:0226" => "ion exchange chromatography",
    "MI:0227" => "reverse phase chromatography",
    "MI:0231" => "mammalian protein protein interaction trap",
    "MI:0232" => "transcriptional complementation assay",
    "MI:0255" => "post transcriptional interference",
    "MI:0369" => "lex-a dimerization assay",
    "MI:0370" => "tox-r dimerization assay",
    "MI:0397" => "two hybrid array",
    "MI:0398" => "two hybrid pooling approach",
    "MI:0399" => "two hybrid fragment pooling approach",
    "MI:0400" => "affinity technology",
    "MI:0401" => "biochemical",
    "MI:0405" => "competition binding",
    "MI:0406" => "deacetylase assay",
    "MI:0410" => "electron tomography",
    "MI:0411" => "enzyme linked immunosorbent assay",
    "MI:0415" => "enzymatic study",
    "MI:0416" => "fluorescence microscopy",
    "MI:0419" => "gtpase assay",
    "MI:0420" => "kinase homogeneous time resolved fluorescence",
    "MI:0423" => "in-gel kinase assay",
    "MI:0424" => "protein kinase assay",
    "MI:0425" => "kinase scintillation proximity assay",
    "MI:0426" => "light microscopy",
    "MI:0428" => "imaging technique",
    "MI:0432" => "one hybrid",
    "MI:0434" => "phosphatase assay",
    "MI:0435" => "protease assay",
    "MI:0437" => "protein three hybrid",
    "MI:0440" => "saturation binding",
    "MI:0508" => "deacetylase radiometric assay",
    "MI:0509" => "phosphatase homogeneous time resolved fluorescence",
    "MI:0510" => "homogeneous time resolved fluorescence",
    "MI:0511" => "protease homogeneous time resolved fluorescence",
    "MI:0512" => "zymography",
    "MI:0513" => "collagen film assay",
    "MI:0514" => "in gel phosphatase assay",
    "MI:0515" => "methyltransferase assay",
    "MI:0516" => "methyltransferase radiometric assay",
    "MI:0655" => "lambda repressor two hybrid",
    "MI:0657" => "systematic evolution of ligands by exponential enrichment",
    "MI:0678" => "antibody array",
    "MI:0695" => "sandwich immunoassay",
    "MI:0696" => "polymerase assay",
    "MI:0726" => "reverse two hybrid",
    "MI:0727" => "lexa b52 complementation",
    "MI:0728" => "gal4 vp16 complementation",
    "MI:0809" => "bimolecular fluorescence complementation",
    "MI:0813" => "proximity enzyme linked immunosorbent assay",
    "MI:0824" => "x-ray powder diffraction",
    "MI:0825" => "x-ray fiber diffraction",
    "MI:0827" => "x-ray tomography",
    "MI:0841" => "phosphotransferase assay",
    "MI:0870" => "demethylase assay",
    "MI:0872" => "atomic force microscopy",
    "MI:0879" => "nucleoside triphosphatase assay",
    "MI:0880" => "atpase assay",
    "MI:0887" => "histone acetylase assay",
    "MI:0889" => "acetylase assay",
    "MI:0892" => "solid phase assay",
    "MI:0894" => "electron diffraction",
    "MI:0895" => "protein kinase A complementation",
    "MI:0899" => "p3 filamentous phage display",
    "MI:0900" => "p8 filamentous phage display",
    "MI:0905" => "amplified luminescent proximity homogeneous assay",
    "MI:0916" => "lexa vp16 complementation",
    "MI:0920" => "ribonuclease assay",
    "MI:0921" => "surface plasmon resonance array",
    "MI:0946" => "ping",
    "MI:0947" => "bead aggregation assay",
    "MI:0949" => "gdp/gtp exchange assay",
    "MI:0953" => "polymerization",
    "MI:0968" => "biosensor",
    "MI:0969" => "bio-layer interferometry",
    "MI:0972" => "phosphopantetheinylase assay",
    "MI:0976" => "total internal reflection fluorescence spectroscopy",
    "MI:0979" => "oxidoreductase assay",
    "MI:0984" => "deaminase assay",
    "MI:0989" => "amidase assay",
    "MI:0990" => "cleavage assay",
    "MI:0991" => "lipid cleavage assay",
    "MI:0992" => "defarnesylase assay",
    "MI:0993" => "degeranylase assay",
    "MI:0994" => "demyristoylase assay",
    "MI:0995" => "depalmitoylase assay",
    "MI:0996" => "deformylase assay",
    "MI:0997" => "ubiquitinase assay",
    "MI:0998" => "deubiquitinase assay",
    "MI:0999" => "formylase assay",
    "MI:1000" => "hydroxylase assay",
    "MI:1001" => "lipidase assay",
    "MI:1002" => "myristoylase assay",
    "MI:1003" => "geranylgeranylase assay",
    "MI:1004" => "palmitoylase assay",
    "MI:1005" => "adp ribosylase assay",
    "MI:1006" => "deglycosylase assay",
    "MI:1007" => "glycosylase assay",
    "MI:1008" => "sumoylase assay",
    "MI:1009" => "desumoylase assay",
    "MI:1010" => "neddylase assay",
    "MI:1011" => "deneddylase assay",
    "MI:1016" => "fluorescence recovery after photobleaching",
    "MI:1019" => "protein phosphatase assay",
    "MI:1024" => "scanning electron microscopy",
    "MI:1026" => "diphtamidase assay",
    "MI:1030" => "excimer fluorescence",
    "MI:1031" => "protein folding/unfolding",
    "MI:1036" => "nucleotide exchange assay",
    "MI:1037" => "Split renilla luciferase complementation",
    "MI:1038" => "silicon nanowire field-effect transistor",
    "MI:1087" => "monoclonal antibody blockade",
    "MI:1088" => "phenotype-based detection assay",
    "MI:1089" => "nuclear translocation assay",
    "MI:1111" => "two hybrid bait or prey pooling approach",
    "MI:1112" => "two hybrid prey pooling approach",
    "MI:1113" => "two hybrid bait and prey pooling approach",
    "MI:1137" => "carboxylation assay",
    "MI:1138" => "decarboxylation assay",
    "MI:1142" => "aminoacylation assay",
    "MI:1145" => "phospholipase assay",
    "MI:1147" => "ampylation assay"

If you also want those that might not be direct binary but simply in the same complex, include:
      "MI:0019" => "coimmunoprecipitation",
      "MI:0006" => "anti bait coimmunoprecipitation",
      "MI:0007" => "anti tag coimmunoprecipitation",
      "MI:0858" => "immunodepleted coimmunoprecipitation",
      "MI:0096" => "pull down",
      "MI:0963" => "interactome parallel affinity capture",
      "MI:0676" => "tandem affinity purification"


Samtools/bcftools/tabix is available in R via Rsamtools.
SEQprocess also provides wrappers for different tools such as GATK and bowtie.
More utility packages are noted here.


Is there a way to quickly verify the presence of some SNPs in Fastq files?

Try this pipeline:

Clean raw fastq files (trimmomatic)
Align clean fastq files to the reference genome (dna-to-dna aligner of your choice)
Convert the alignment sam file to bam (samtools)
Call SNPs from bam file (samtools)

Here's a start tutorial on this: http://ged.msu.edu/angus/tutorials-2013/snp_tutorial.html


NCBI SRA database sample: control vs test

You can download the information using this link:

I did it in R:

tab = read.delim("https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=SRA059451&result=read_run&fields=run_accession,experiment_title,sra_ftp,sample_alias")

head(tab)
  run_accession
1     SRR585697
2     SRR585700
3     SRR585726
4     SRR585727
5     SRR585728
6     SRR585729
                                                                experiment_title
1 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 171
2 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 238
3 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 241
4 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 242
5 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 243
6 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 245
                                      sra_ftp
1 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585697
2 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585700
3 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585726
4 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585727
5 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585728
6 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585729
                                    sample_alias
1     Symptomatic atherosclerosis gut metagenome
2 Symptomatic atherosclerosis gut metagenome 238
3 Symptomatic atherosclerosis gut metagenome 241
4 Symptomatic atherosclerosis gut metagenome 242
5 Symptomatic atherosclerosis gut metagenome 243
6 Symptomatic atherosclerosis gut metagenome 245

Normally the information will be here, but I cannot find what does metagenome 238 etc refers to. No table in the paper that can link sample to phenotype..

Most likely you need to contact the authors. This list might be useful though.


Pathway enrichment analysis with or without custom background

I would defninitely go with a custom background which should be all the proteins that were actually measured in your experiment since you should only take into account what was measured.


extracting last few lines of a growing bgzip compressed .vcf.gz file

27GB is not huge. You can zcat the file and see the last few lines. It will take ~10 minutes. It is theoretically possible to peek a bgzf file, but no such tools are available.

By the way,

  I am running bcftools merge to create a single .vcf.gz file out of a collection of chromosomal .vcf.gz files from the 1000g project.

Probably you should use bcftools concat, not bcftools merge. merge combines VCFs from different samples. concat concatenates VCFs for individual chromosomes; each input VCF contains the same set of samples.


How to design DESeq2 LRT model with individuals nested in 2 levels?

From what I can gather, you want to account for the effect of individual, nested within region. That is, you want to see after accounting for these, is there a consistent effect for Injection:Social across all conditions.
So you set up the model like this:
m1 <- model.matrix(~ ind.n*Region + Injection + Social + Injection:Social,data=..)

The last term should be Injection:Region and you can just use the waldTest (default) in DESeq2 for this term.
What does the terms do? ind.n*Region is the equivalent of ind.n + Region + ind.n:Region , and with this you effectively get an effect for every region in every individual.
Why don't we need the Injection:ind.n or Social:ind.n or Social:Region. These terms indicate the effect of Injection or Social can vary by individuals or regions. Most likely introducing too many parameters when you are interested in a common effect. Also you do not have the replicates or samples to distinguish this effect from region or other effects.
Since you provided an example we can run DESeq2 and you can see how the results look like:
mat = counts(makeExampleDESeqDataSet(n=1000,m=48))
dds = DESeqDataSetFromMatrix(mat,df,~ ind.n*Region + Injection + Social + Injection:Social)

dds = DESeq(dds)

resultsNames(dds)
 [1] "Intercept"            "ind.n_S2_vs_S1"       "Region_NS_vs_HY"
 [4] "Region_TN_vs_HY"      "Injection_S_vs_L"     "Social_KF_vs_ISO"
 [7] "Social_KFC_vs_ISO"    "Social_NF_vs_ISO"     "ind.nS2.RegionNS"
[10] "ind.nS2.RegionTN"     "InjectionS.SocialKF"  "InjectionS.SocialKFC"
[13] "InjectionS.SocialNF"

The terms you need are "InjectionS.SocialKF","InjectionS.SocialKFC", "InjectionS.SocialNF", and you can look at each of them:
head(results(dds,name="InjectionS.SocialNF"))
log2 fold change (MLE): InjectionS.SocialNF
Wald test p-value: InjectionS.SocialNF
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1  9.9811166787259   1.25304112986447 0.819806376919295   1.52845984752303
gene2 30.3449455820337 0.0329442893152027 0.705199688255367 0.0467162562092241
gene3 3.83223545055379    1.0281136369045  1.64095596190233  0.626533350543196
gene4  11.232305747171  0.595738624408923   0.8243883031544  0.722643227868976
gene5 6.70950627004097  0.756449993378065   1.0631622863378  0.711509430967263
gene6 26.1431134888287 -0.854784518963918 0.625714541243558  -1.36609342219393
                 pvalue              padj
              <numeric>         <numeric>
gene1 0.126398405431826 0.978671002658464
gene2 0.962739373909937 0.999897888026606
gene3 0.530965168963838 0.978671002658464
gene4 0.469899103018657 0.978671002658464
gene5 0.476768608734069 0.978671002658464
gene6 0.171909642630577 0.978671002658464

As mentioned, you can do a LRT if you want to test all the Injection:Social term interactions terms at one go, that is, the null hypothesis is that all of them are zero:
dds = nbinomLRT(dds,reduced=~ ind.n*Region + Injection + Social)
results(dds)

Usually the individual terms are more intuitively than testing all of them are zero, but you might have a special need for this.


You can have a look at cutadapt. It is capable of quality trimming for both ends as you can read here.


The following should work:

cat(file="out.fa", paste(paste0(">",names(dnabin_to_fasta)),
          sapply(dnabin_to_fasta, paste, collapse=""), sep="\n"), sep="\n");

A bit of explanation, starting from the middle outwards:

Collapse the letter arrays inside dnabin_to_fasta into single strings
Append ">" to the front of the names of each sequence
Concatenate the results of steps 1 & 2 together, separating with a line break
Write each concatenated string to a file, separating by a line break


Common QC parameters and threshold in human WGS pipeline

I would suggest you look at the website for the excellent MultiQC tool (https://multiqc.info/).  We use it in a clinical setting to collate all the QC statistics from our pipelines into a single report (It displays results from multiple QC programs such as FastQC, Picard, bcl2fastq, VerifyBAMid etc).  The homepage has a selection of example reports - clicking on "Whole Genome Seq" shows a report with the metrics most labs would expect from their WGS pipeline.

If you dig around in their GitHub repo (https://multiqc.info/) you will get a good idea of the default thresholds set for flagging results - there should be a config file somewhere specifying all these values.


Input from multiple channel, merge and use in nextflow

I don't think you can unpack tuples using the each qualifier like that. You can, however, unpack them in the script block:
process test {

    input:
    each triple from para1

    script:
    def (m, k, l) = triple

    """
    echo "${m} ${k} ${l}"
    """
}

Consider also that the merge operator is deprecated and will be removed from future Nextflow releases. Instead, you may need to join by index, for example:
def indexedChannel( items ) {
    return Channel.from( items.withIndex() ).map { item, idx -> tuple( idx, item ) }
}

k_ch = indexedChannel( [ 15, 20, 21 ] )
m_ch = indexedChannel( [ 'a', 'b', 'c' ] )
l_ch = indexedChannel( [ 1, 2, 3 ] )

m_ch
    .join( k_ch )
    .join( l_ch )
    .map { idx, m, k, l -> tuple( m, k, l ) }
    .into { para1; para2 }


Usually it is more problematic to work with older versions of data with newer protocols than the other way around.

In this case, I think that the features annotated in the (newer) bed file should be more accurate than the older one. Unless you specifically want to compare your experiment with a previous experiment performed with a bed file from 2012 you'll be fine using the newer one.

However make sure that you are correctly using the right file and that you report (or record) each version of protocol and data used in case some doubts arise.


Can you query the RCSB PDB API with a ligand and return the IDs of all entries in which it appears as a free ligand?

I'm not sure about the RCSB, but you can get such a list from the PDBe REST API, as documented in the compounds section:

http://www.ebi.ac.uk/pdbe/api/pdb/compound/in_pdb/ADN

Alternatively, from an older service:

http://www.ebi.ac.uk/pdbe-srv/pdbechem/PDBEntry/download/ADN


You can easily do it via Notepad++ or even Excel. Just don't forget to save the output file as a .csv file (comma-separated value).
The info inside should look like (my working example):

id,state
lar_control141,control
lar_control142,control
lar_control143,control
lar_probiotic142,probiotic
lar_probiotic143,probiotic
lar_probiotic144,probiotic
lar_probiotic145,probiotic
lar_probiotic146,probiotic

That should work


In my opinion, it is not very reliable. LiftOver is very limited in terms of transformations it can support. The LiftOver Chain format can capture only matching regions in the same order. It means that it can account for indels, but even simple structural variations become problematic.

For instance, when a newer assembly is available, it is usually a recommended practice to remap all the reads rather than transform the existing alignments.


Getting structure and sequence related to PDB IDs

You can download seqs and structures based on a list of PDB ids using http://www.rcsb.org/pdb/download/download.do#FASTA


How can i get data set of illness and not illness people?

The data set will be available on Genbank, you simply search for the "authors surname and HBV" and their entire data set should be revealed.

The HBV genome is 5Kb, the smallest viral genome so its not alot of data.

I am confused by "Genotype B nucleotide sequence and non-HBV nucleotide sequences". What I suspect is the investigators looked at HBV genotype B in comparison to HBV genotypes A, C to I. Essentially they wanted to assess the relationship between genetic variation and pathogenesis, which is an attempt to assess virus-induced cancer. Obviously I don't know the paper but I'm aware of the basic research front in this area. HBV has a very wierd evolution and I personally wouldn't want to investigate it.


How to download RNAseq gene expression data from GTEx

Gene read counts are the number of reads or fragments that are mapped to each gene. So it is an integer and the values are not directly comparable between samples due to differences in sequencing depth. TPM stands transcripts per million. It is the normalized gene expression level. Basically, you first normalize read counts by gene length and then normalize the library size to 1e6 in all libraries. As for which one should you use, it depends on what analysis you want to do and what tool you are going to use.
I have not noticed that there are duplicated genes in those files, can you show an example?
I do not think you can download expression from one tissue, but you can always subset. The whole dataset is not very large.


merge the matrix from computematrix of deeptools

You're trying to cbind two matrices of different sizes, so computeMatrix is throwing an error since it's unclear which row goes with which. You need to ensure that matrices have exactly the same number of rows and that the Nth row in each match (i.e., they should be bound together with cbind).


what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?

The best set I'm aware of can be downloaded from https://gnomad.broadinstitute.org .


Tools to do VCF to MAF and MAF to VCF conversion?


Error processing run file with bcl2fastq - Unable to open '.../RunInfo.xml' file for reading

I repro the error.
ls /ml

ls: cannot access '/ml': No such file or directory

If run
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

I suspect it will work.
At present you are attempting to call a directory /data, which doesn't exit. If you type ls /etc that will work. The difference between ./data and /data in a file structure is alot, but on a command line it appears minimal.

Wait ... if that doesn't work .. because I suspect the calling directory is data then the following will,
bcl2fastq --runfolder-dir ./201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

OR
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./data/SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

It depends on the path where you are making the call. You could sort this out using $PWD or just echo $PWD then copy and paste the path name.


At what stage of a transcriptome assembly is it better to perform read contaminant filter?

First, to answer your question about mapping to a low-quality reference:

1. Mapping

For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome assembly does not contain any scaffolds from contaminating species like bacteria, algae, et cetera.

2. Reasons to not BLAST your raw reads

There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in #4, then reassembling, may work better.

3. Using multiple individuals of highly heterozygous species in assembly.

Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.

4. A strategy to filter transcripts based on mapping and BLAST

Given #1, #2, and #3 above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.

I would try this:

Assemble the transcriptome of each individual separately.

This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.

Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.

Use somthing like cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta.

Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.

Use samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam to only keep the transcripts that did not mapped to the reference. -f 4 is the magic flag here that keeps things that didn't map.

Make a list of all the transcripts that did not map to the reference...

samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt

...and get only the transcripts from your fasta file that are potential contaminants.

seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta

BLAST all of the unmapped.fasta and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this known_contaminants.fasta
Map all reads from all individuals to the known_contaminants.fasta.
Get all the reads that map to known contaminating transcripts, known_contaminants.fasta using:

samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt

Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
Reassemble the filtered fastq files.

This strategy has the following benefits:

Does not require weeks/months of compute time to BLAST 2 billion reads.
Ensures that all reads from the reference genome are left untouched.
Ensures that all reads from known contaminants are removed.
Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
Leverages contaminants from different individuals to filter reads from all individuals.


Genotyping with Illumina HumanOmni1-Quad before whole genome sequencing

The Illumina HumanOmni1-Quad beadchip is a microarray device consisting of 1,140,419 markers which have derived from the 1,000 genomes project. The markers chosen are is high-value regions of the genome: ABO blood typing SNPs, cSNPs, disease-associated SNPs, eSNPs, SNPs in mRNA splice sites, ADME genes, AIMs, HLA complexes, indels, introns, MHC regions, miRNA binding sites, mitochondrial DNA, PAR, promoter regions, and Y-chromosome.

The markers come together to form a sort of individual fingerprint. Germline SNPs are heritable and so can be traced up the line. Some SNPs are more prominent in subpopulations, something that Ancestry.com and other similar services leverage to predict ethnicity.

This paper utilized genomic DNA extracted from the Singapore BioBank. BioBanks tend to record age, lifestyle, family history, and health information alongside the genomic DNA. I (conservatively) assume that the Singapore BioBank at least has familial relations and gender information tied to the DNA. My interpretation is that the ivestigators SNP genotyped using the Illumina Beadchip the genomic DNA from the bank and compared their results with the "ground truth" provided by the genetic bank.

Illumina brochure for the HumanOmni1-Quad beadchip.


Why is a PacBio read length larger than the aligned reference region?

Q1
Probably not - this is probably just displaying the contents of the input (I presume a bam file) exactly as they are encoded.
Q2
The splicing looks a little weird, but the exon-mapping pattern in the 5' part of the gene makes it seem like RNA rather than DNA. I think that the large intron is what makes you think that it is DNA, correct?

Explanations for weird splicing

The weird read is a rare splice variant but a true mRNA.
The weird read is not completely matured and has an intron that has not yet been spliced. The IsoSeq library prep protocol picked this up somehow and it ended up being sequenced.

Q3

Explanations for why Seq length: 6751 and Seq length on ref: 5267

There is probably a large amount of clipping on the 5' or 3' end of the read. This could be polyA or some other transcribed-but-not-yet-completely-matured seq. Or, the extra sequence could be sequencing error, PacBio adapters, or a chimeric molecule.

things to try
How to check what is going on:
I'm guessing that you had to upload a bam file to JBrowse to make the IsoSeq track. Use the sequence name that you got from Jbrowse and run the following command to extract the alignment info from the bam file:
samtools view /path/to/your/isoseq/file.bam | grep 'weird_read_id'

Look at the cigar string to see what it says about 5' or 3' clipping and corroborate with the sequence data.


conda doesn't install latest version of snakemake

Python (whether installed with conda or not) will prefer packages installed in your home directory unless you tell it not to. This is a little known secret about using a virtualenv, actually, since internally those tell python to ignore everything in your home directory. Within conda, there's long been a debate about whether python (and also R) should continue looking for packages in the normal places or only within a conda environment. Since it would likely break many people's workflows, it's been decided to allow the normal language features of looking for packages not just within conda envs but also in the normal system and user locations. For python in particular you can disable this as follows:

PYTHONNOUSERSITE=True snakemake ...

Then only the packages installed within your conda environment will be used.