Aligning ChIP-Seq reads to repeats for downstream peak analysis

This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips.

  However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

If this dataset is also PE and the alignment is outputting either uniquely mapped reads or a single (random) position for multimappers, it is pretty much as good as it gets. A lot of repeat regions are not exact copies, so using PE there is a good change that one read of the pair mapping to a single position is helping with mapping of the other read.

  I need to know reasonably confidently where precisely those repeats are located

It is probably a matter of phrasing, but what one needs is to know the precise genomic location of the reads. The location of the repeat elements can be obtained from the Repeat Masker track of UCSC, but you probably know this already.

  With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I think bowtie already does this. Did you map the read pairs together or separately?

Note that:

There will always be a number of multimapping reads (fragments) that will be discarded. The information of how many are being discarded for your samples would be helpful.
Mapping can be improved using both reads and tunning bowtie settings (see bellow)

  Finding a valid paired-end alignment where both mates align to
  repetitive regions of the reference can be very time-consuming. By
  default, Bowtie avoids much of this cost by imposing a limit on the
  number of "tries" it makes to match an alignment for one mate with a
  nearby alignment for the other. The default limit is 100. This causes
  bowtie to miss some valid paired-end alignments where both mates lie
  in repetitive regions, but the user may use the --pairtries or -y
  options to increase Bowtie's sensitivity as desired. source

  The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced.

As indicated above, the aligner will take care of that sort of. It will gave the mapping information of each read pair, which is then used by MACS.

Other tools

I tested ShortStack which uses the coverage of uniquely mapped reads to place multimapping reads for ChIP-seq, and the results were mixed. That said, it was SE data and the code has been through some changes, so it be worth a try.

When requested, the authors of RepEnrich provided me with the scripts for an utility called RepConsensus which maps reads to a consensus reference of repeat regions. The caveats are:

mappings are better for elements whose sequences have not diverged a lot yet
it will tell the position of the protein in the repeat sequence, but no information about the surrounding genomic regions because there isn't any.

Whatever you do, do not use all alignments, that is double counting of reads, and think very carefully before using custom methods. A blog post detailing some of the potential issues can be found here (I wrote it).


What is the best argument for maximum parsimony method in phylogenetic tree construction?

I'd say maximum parsimony approaches assume genetic similarity is unlikely to occur by convergence, but likely a result of divergence. If this assumptions holds, finding the tree configuration that minimizes the number of unlikely events makes sense.


Highly variable genes analysis by DESeq2

The first method (resSig$log2FoldChange) is for getting the most differentially expressed genes. These are the genes with the biggest differences between specific groups you pre-define.

The second method (rowVars(assay(rld))) is for getting the highly variable genes. These are the genes most variable across all samples regardless of which samples they are.


Aligning sequence and comparing it against primer

If I understand you correctly, after your most recent edit, you have 5,000 sequences that are 63 bp long and in each you simply want to replace the nucleotides from 21 through 43 with - unless that nucleotide does not match the primer. Since - is normally used to represent gaps, I recommend using . instead for matching nucleotides.

I know you're looking for an R solution, but I don't usually use R for string manipulations and don't have time at the moment to figure it out. So, assuming this might be urgent for you, below is a Python solution. This uses your example data which shows 11 nucleotides before the primer and 9 after it (so you'll need to adjust it to work with your actual case in which 20 nucleotides are before and after the primer). In Python, strings are arrays of characters, by the way.

Here are your example sequences are in the file seqs.txt:

ATGACCATACGATGTCGAAGCAAATATCAATTGAATCGTAGGG
ATGAGCATACCATGTGAAAGCTAATATCAATTGATTCGAACGC

Here is the code that gives you the expected output:

#!/bin/env python
import sys
primer = 'ATGTGGAAGCAAATATCAAATGA'

print('-' * 11 + primer + '-' * 9 )
with open('seqs.txt','r') as f:
    for gene in f:
        sys.stdout.write(gene[0:11])
        for gene_index in range(11,34):
            primer_index = gene_index - 11
            if gene[gene_index] == primer[primer_index]:
                sys.stdout.write('.')
            else:
                sys.stdout.write(gene[gene_index])
        print(gene[34:])

Running this gives you:

-----------ATGTGGAAGCAAATATCAAATGA---------
ATGACCATACG....C..............T...ATCGTAGGG
ATGAGCATACC.....A....T........T...TTCGAACGC


How to export peptide-shaker data from commandline?

You have to start PeptideShaker calling java with the -cp argument an input project file and provide a comma/space-separated list of integers for the reports you want to export:

java -cp PeptideShaker-1.16.37.jar eu.isas.peptideshaker.cmd.ReportCLI -in input.cpsx -reports "0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10"`

Here is the list of reports you can export:

-reports                   Comma separated list of types of report to export.
                           0: Certificate of Analysis,
                           1: Default Hierarchical Report,
                           2: Default PSM Phosphorylation Report,
                           3: Default PSM Report,
                           4: Default PSM Report with non-validated matches,
                           5: Default Peptide Phosphorylation Report,
                           6: Default Peptide Report,
                           7: Default Peptide Report with non-validated matches,
                           8: Default Protein Phosphorylation Report,
                           9: Default Protein Report,
                           10: Default Protein Report with non-validated matches,
                           11: Extended PSM Report,
                           12-n: Your own custom reports.

I don't know whether you can also use the SearchGUI output zip to generate this files. The program accepts zip files as input, but I think that means PeptideShaker project that are exported to zip files, but not SearchGUI output zips. So, in the end you might need to run the GUI at least one time, unfortunately.


How do I generate a color-coded tanglegram?

I think you can try dendextend, in this manual there is an example of coloring the branches. I don't think it is exactly like your coloring, but with a little tweaking you might get your colorscheme in there.

The manual mentions an argument called color_lines for the function  tanglegram():

# The `which` parameter allows us to pick the elements in the list to compare
iris_dendlist %>% dendlist(which = c(1,4)) %>% ladderize %>%
   # untangle(method = "step1side", k_seq = 3:20) %>%
   set("rank_branches") %>%
   tanglegram(common_subtrees_color_branches = TRUE)

According to the manual, the code above produces an image like this:

Try to see if you can fit your color palette in there.


Clarification regarding (deepTools) bamCoverage output

Yes
The total RPKM for the region in the bedgraph file is the value already present in the bedgraph file. If you want to estimate the RPKM for a larger region composed of multiple bins then you can average the per-bin RPKM values. Note that adjacent bins with the same value are merged by deepTools, so you will need to account for this.

Note that RPKM is NOT a robust normalization and that I should remove it from deepTools so people stop using it.


GSEA preranked different results opposing phenotypes

What do you mean by having the same results but opposite direction?
What GSEA does is that it goes through the gene list from the top to bottom, whenever it encounters a gene that belongs to a gene set, it adds a positive score to that gene set, and if not, it penalizes the score. That is what you see in the green line in the results.
Now imagine that you reverse the order of your input gene list, it is obvious the results are going to be different if the program still goes through the list from top to bottom.
I think as a general rule, whenever you change the order of input gene list, results are going to be different.


Assuming the lists are named, and you want to fetch the element with the same name from each list, I'm guessing something like this would work (or at least get you closer to a solution):
v <- lapply(names(x), function(var){
   forest(x[[var]], slab=(my_list[[var]])$StudyID)
})

When creating images that don't have a multi-page format (i.e. anything except PDF), and not changing the image name with each save, you need to add in a format specifier to tell R what to do with new created images. For example, this will add a three digit zero-padded number for each new TIFF image:
tiff(filename = "~/Desktop/image_%03d.tiff", width = 30, height = 30, units = "cm", res = 1200)


Is there a standard way to clean a PDB file and re-number its residues?

Several years ago we have written a tool for this: http://kodomo.fbb.msu.ru/~avkitex/projects/pdbparser.html
https://github.com/avkitex/pdbparserTools

The description is only in Russian but you can contact the author.


Total read count of shotgun affected by few highly abundant sequences

That sentence in the Anders and Huber paper is a bit misleadingly worded (I don't think that was their intent). What they mean is that highly expressed and highly differentially expressed genes can out-compete the majority of the other genes, which thereby leads to them appearing to have lower expression if one uses total read counts for normalization (basically, what you said in the second paragraph).

An alternate interpretation of what they wrote (likely the one you understandable understood to have been their meaning) is that some highly expressed genes will inflate the total number of reads a sample actually gets. I have yet to see this in practice, but it certainly COULD happen. There are technical biases in sequencing such that some sequences (often due to length and GC content) are easier to sequence on common machines. A good example of this would be PhiX, which forms clusters much more easily on Illumina machines than real samples. I expect that this is relatively rare in practice, but using a more robust normalization would still help in protecting against it.


How to count RNA sequencereads using custom made python scripts?

What software are you using for the next step?  DESeq in R will do this filtering for you in one line.
I think a more typical way to do this would be to make your entire gene count matrix, and then remove rows where the max value for that row is < 5.  Your way, you have to go through each file twice; once to know which genes to keep, and then again to actually remove lines.


Not clear what we are trying to do, but here is the solution using lapply that should replicate your existing forloop.

# example input data
snp_ids <- read.table(text = "
snp_id  chr  pos
Chr01__912   1  912
Chr01__944   1  944
Chr01__1107  1 1107
Chr01__1118  1 1118
Chr01__1146  1 1146
Chr01__1160  1 1160", header = TRUE)

# window sizes to loop through
win_size <- c(15000, 30000, 50000, 100000)

res <- cbind(snp_ids,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 as.numeric(paste(snp_ids$chr, ceiling(snp_ids$pos/w), sep = "."))
               )))

# result
res
#        snp_id chr  pos window_15000 window_30000 window_50000 window_1e.05
# 1  Chr01__912   1  912          1.1          1.1          1.1          1.1
# 2  Chr01__944   1  944          1.1          1.1          1.1          1.1
# 3 Chr01__1107   1 1107          1.1          1.1          1.1          1.1
# 4 Chr01__1118   1 1118          1.1          1.1          1.1          1.1
# 5 Chr01__1146   1 1146          1.1          1.1          1.1          1.1
# 6 Chr01__1160   1 1160          1.1          1.1          1.1          1.1

Edit: If we really want to use forloop, then try below:

for (i in win_size){
  windows <- ceiling(snp_ids$pos/i)
  snp_ids[, paste0("window_", i)] <- as.numeric(paste(snp_ids$chr, windows, sep = "."))
}


You can use pymolrc and create aliases to your scripts.

If your script is located e.g. in

/Users/tillyboy/pymolscripts/myscript1.pml

Just enter this line to your pymolrc:

alias myscript1, run /Users/tillyboy/pymolscripts/script1.pml

Now myscript1 is available for all your Pymol sessions independent of your current working directory.


probability of an entire sequence, not just one value

The author writes in the quote that $O$ is to be taken for the whole set of positions $O_{i}$.
In the extracts that you provide there is no reason to believe that it is only for a sequence of length 1. The author writes "the observed sequence $O = O_1, O_2...$"
Thus as you say $Pr(sequence) = Pr(O) = Pr(O_1) * Pr(O_2)...$ according to my understanding of the quotes.


You are using the argument values incorrectly. Argument values are sample attributes, not full character paths. Hence, "IU" is a sample attribute, "gencode_mouse_m13" is a sample attribute, so why would you think "/scratch..." would be treated any differently?

The pipeline interface maps sample attributes to command-line arguments. It is modulated by sample (different values per sample). What you're trying to do is pass a constant string to an argument.

For this, there is a pipeline_args section of your project config. Please read the docs on project config.


Converting hmmer --tblout to pandas dataframe

The manual says fields in a tabular output file are space-delimited. It should be fairly straightforward to load a space-delimited file in pandas.

From the manual:

  --tblout <f> Save a simple tabular (space-delimited) file summarizing the per-target output, with one data line per homologous target model found.


How to interpret duplication from MultiQC/FastQC?

I would agree with your statement.  Typically the ratio of unique to duplicated reads is a function of input library complexity (number of input molecules) verse sequencing depth (the amount of molecules sequenced).

In the example you gave you have about a 2:1 ratio of duplicates to unique reads which suggests that these libraries have been sequenced to their fullest and additional sequencing would yield few new unique reads.


You just need to check the flag for suplementary mapping. If there is suplementary mapping then use the coordinates; if they are all the same the seq is the same and you just have one event, if different reads have different suplementary mapping locations then there's multiple events. The problem is more complex than this though so just use an existing tool.


How does one convert a log-likelihood substitution matrix into a probability matrix?

UPDATE2
@chubs wants to compare different clades. This is possible within recent releases of ConSurf:

It is additionally quite possible to use ConSurf's default outputs to get per-position conservation scores that can be used for post-processing in any fashion, as well as prefab scripts that can be used with PyMol to make custom visualization:

For all cases, ConSurf creates the following outputs:

The sequence and MSA colored by ConSurf conservation scores.

A text file that summarizes for each position the normalized score calculated, the assigned grade (1-through-9), the reliability estimation (for the Bayesian method) and the amino acids/nucleotides observed in the respective MSA column.

The sequences selected for the MSA and the MSA constructed (unless those files were uploaded by the user).

A file with the frequency of each amino acid/nucleotide observed in each column of the MSA.

The evolutionary tree, which was calculated by the server or uploaded by the user, together with the MSA are shown using the WASABI platform (Veidenberg et al., 2015).  Moreover, using WASABI, a user can select a subtree containing a fraction of the homologous sequences and conduct a follow-up ConSurf analysis with these selected sequences. To refine a ConSurf analysis to a selected subtree, a user can choose any internal node on the phylogenetic tree and open a WASABI menu by using a right mouse click. Selecting the option 'run ConSurf on subtree' will issue a new window with the new ConSurf run for the selected subtree sequences (see example in Figure 3).

UPDATE
It seems that what @chubs really wants to do is map protein conservation onto a structure. There are various tools for doing this (here, here). ConSurf is the older and more established tool, using evolutionary conservation to rule out trivial apparent associations due to phylogenetic non-independence.
here is the  somewhat related 3DPatch workflow, which uses profile HMM information as a measure of conservation:

Here is the output that they get:

Likewise ConSurf gets similar looking outputs:

ORIGINAL ANSWER
I think that more information is needed here to answer this question properly (what is your data? what is your ultimate goal? what are the things you've already tried?), but we can say a few things right away.
I would strongly suggest looking into substitution models, which are literally substitution probability matrices. Note that these are probabilities per unit time, which is to say that they are calibrated to e.g. evolutionary rates along a specific phylogenetic tree. I would therefore strongly suggest incorporating a phylogenetic model of some kind.
If you have e.g. a multiple sequence alignment, you can estimate a phylogenetic tree by one of the many tree estimation methods, which will give you a substitution probability matrix "for free" as this estimation is taken care of in the process of tree estimation. These will again be substitution probabilities per unit time; you will honestly not be able to get rid of the time-dependence here, as that's how evolution works.
I am a little behind the times, but some tools you might look into for tree estimation are PhyML and FastTree. If you only have sequences and haven't aligned them yet, you can use a tool like MUSCLE, ClustalO, or MAFFT.
BLOSUM and related models are derived from specific slices of protein sequence identity (evolutionary history) calibrated to certain scopes, it's unlikely that they are appropriate to your case, unless your data fits inside those rather narrow slices.


How do I validate a single sample ArrayCGH result?

Why not just good old qPCR? That's (A) quick, (B) cheap and (C) easy to analyze. If you care about the exact location of the break point (I'm guessing from the context that you don't), then targeted sequencing with a custom capture kit would work.

Regarding validating the biological relevance of this, there are multiple (parallel) routes one can take. Firstly, screen unaffected family members for this same alteration. If you find this in unaffected individuals then it's obviously not the causative alteration. Ideally, one would also do either a cell-line experiment or a mouse (or other model organism) experiment to see if you can reconstitute at least some component of the clinical phenotype. This may not always be possible, of course.


How can I calculate coverage at single bases using a bam file?

Note: not yet tested, so there may be some additional fiddling with command line options needed

The per-base depth can be obtained from samtools depth (-a includes zero-coverage positions):

samtools depth -a in1.bam > depth_in1_both.tsv

To split this by forward and reverse, you can use an initial pipe through samtools view to exclude or include reverse-complement mappings:

samtools view -b -F 0x10 in1.bam | samtools depth -a > depth_in1_fwd.tsv
samtools view -b -f 0x10 in1.bam | samtools depth -a > depth_in1_rev.tsv

To filter for particular regions, this can be pre-piped through bedtools intersect:

bedtools intersect -a in1.bam -b regions.bed | \
  samtools view -b -F 0x10 in1.bam | \
  samtools depth -a > depth_in1_regions_fwd.tsv

Combining both, forward, and reverse is a bit trickier. I'd recommend doing it in a more manual fashion in R. But if you want to keep to the command line...

If you are absolutely sure that the 'depth' output will appear in the same order in all files (it might be worth using samtools depth -aa), then the lines can be pasted, with anonymous pipes, with unnecessary fields filtered out using cut:

paste <(samtools depth -aa in1.bam) \
      <(samtools view -b -F 0x10 in1.bam | samtools depth -aa) \
      <(samtools view -b -f 0x10 in1.bam | samtools depth -aa) | \
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

Or equivalently, if using intermediate files that have been generated already:

paste depth_in1_both.tsv depth_in1_fwd.tsv depth_in1_rev.tsv |
  cut -f 1,2,3,6,9 > depth_in1_bfr.tsv

It would also be possible to use join to combine the files (in cases where file output is not consistent), but the first two fields of each file would then need to be combined, because 'join' only joins on the first field... but if you need to get that complex, it'd be better to use something other than pipes on the command line.


Picard CollectGcBiasMetrics ignoring certain chromosomes/sequences

The proper solution is to use a different tool (some of this you could do with computeGCBias in deepTools). But since you don't want to do that, you'll have to manually remove unwanted chromosomes from the BAM file(s):

$ cat foo.awk
BEGIN{split(excludes, excludeList, " ")}
{exclude=0
if($1 == "@SQ") {
    for(ex in excludeList) {
        if(ex == substr($2, 4)) {
            exclude=1
            break
        }
    }
} else { \
    for(ex in excludeList) {
        if(ex == $3) {
            exclude=1
            break
        }
    }
}
if(!exclude) print}

$ excludes="chr1 chr2"
$ samtools view -h alignments.bam | awk -v excludes="$excludes" -f foo.awk | samtools -bo alignments.filtered.bam

The awk script is stripping unwanted chromosomes from the header and the body. You can then run the result through picard.


From my memory of what a statistician told me, a PCA aims to determine independent linear combinations of variables (i.e. genotypes) that account for the most variation in the dataset. With 10 million SNPs, the vectors describing the genotypes of all individuals reside somewhere inside a 10-million-dimension hypersphere, and the first principle component determines the best vantage point / angle from outside that sphere that will create the greatest spread of data when data is projected onto a perpendicular / orthogonal axis. The second principle component is the perpendicular / orthogonal axis that will create the second-greatest spread of data.

Aside: There's some wonderful linear algebra (which I've forgotten) that makes it very quick (computationally) to do this angle calculation, exploiting the nature of linear transformations and orthogonality.

The point of the PCA for GWAS is to identify (and preferably exclude) differences that might be due to the structure of populations, rather than the disease of interest. It's fairly common for researchers to do a PCA on the entire set of genotype data, then use the first 5, 10, or 20 principal components as covariates in the association model.

There's an underlying assumption that the majority of genetic variation is due to population structure - if this is not the case, then the disease would have to be extremely omnigenetic, with almost the entirety of the genome contributing to the disease. The approach also won't deal with all the population structure (because there may be non-linear contributions of genotypes to population structure), but it will get rid of the largest problems.

Additional to this, I'd recommend the following things:

Structure Visualisation

Carry out a STRUCTURE analysis or PCA on the combined cases and controls, and display the results of the first two principle components, labelling, colouring, or otherwise indicating both groups. As an additional check, it may be useful to do addititional PCAs for cases and controls separately.

If they appear different, it suggests that the cases and controls were drawn from separate population clusters, which should be commented on in the discussion of the paper. There's no need to panic and find more data, just make the observation.

Subsampling

Carry out association tests on an equal-sized subsample of both datasets. This will mean that area-under-the-curve statistics can be properly interpreted. If you've got unbalanced datasets, use a little bit less than the smaller of the two groups (e.g. for 5k diseased and 95k healthy, use a random subsample of 4.5k from both groups).

Multiple testing / Internal validation

Run the association tests multiple times (with different subsamples), recording the individuals that were subsampled for future reference / reproducibility.

If an association is only present in one of the many different subsamples, it should be ignored as a spurious result. It means that there is some association that is not generalisable to all (or most) of the cases and controls. These variants might be interesting in themselves on an individual level, but not in a general population.

I explain some of this in more detail in the following preprint (which I'm not planning to submit to any journal, due to it being a bit stale and needing a lot of work to fix up):

Bootstrap Distillation: Non-parametric Internal Validation of GWAS Results by Subgroup Resampling


bedtools intersect on very large .bam file - -sorted confusion

It sounds like bedtools is behaving properly. The bam file is sorted by read name only if you used the -n option in samtools sort. Otherwise, if you used samtools sort without the name-sort option, it is 'ordered' by reference name (which can be arbitrarily ordered, as I'll explain below) and then by position (which is numerically ordered).

I think you are getting confused about 'ordering' the reference names. There are multiple ways to sort names; Some people might want numeric sorting (chr1, chr2, chr3...) and others might want lexicographic sorting (chr1, chr10, chr11...). Because of this, many tools allow you to specify the way you want to order the reference names.

There isn't necessary only one correct way to sort the reference names; what matters is that the order is the same between the files you're trying to intersect. This is exactly the purpose of bedtools intersect -g.

If you want to check the order of reference names in your bam, you can easily do so by checking the header via samtools view -h input.bam. It will say something like this:

@HD     VN:1.3  SO:coordinate
@SQ     SN:chr1     LN:248956422
@SQ     SN:chr10    LN:133797422
@SQ     SN:chr11    LN:135086622
@SQ     SN:chr12    LN:133275309

The @SQ lines tell you how your bam is sorted. My suspicion is that this order is different from the Linux sort or bedtools sort default, which is why your bedtools intersect was failing.


Exon-exon junctions: compare experimental transcripts to reference annotation

You can use the script that comes with minimap2:

# install the k8 javascript shell
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8    # assuming $HOME/bin in your $PATH

# acquire the latest minimap2
git clone https://github.com/lh3/minimap2.git

# print novel/erroneous junctions
minimap2/misc/paftools.js junceval -e anno.gtf rna-seq.sam > output.txt

Some example output:

P       read1    4       chr2    227332588       227340301       [(227332588,227340291)]
P       read2    1       chr9    19376379        19376489        [(19376388,19376493)]
P       read2    3       chr9    19378506        19378707        [(19378514,19378707)]
N       read3    1       chr1    160647335       160647472
N       read3    2       chr1    160647581       160659672

Here column 4–6 gives junctions not in the annotation. Column 7, if not empty, gives the list of annotated junctions/introns that overlap the junction in read alignment. Note that a big fraction of outputted records are errors in read alignment.


Where can I find a list of SNP rs IDs that belong to the X chromosome?

You just searched for X[Chromosome], you didn't specify a species. Presumably, your dataset comes from a specific species, so you should limit your search to that species. For instance, for human, you can search dbSNP for X[Chromosome] human[organism]. Then, click on the send to link at the top left and download the data:

This is veeeeeery slow, however. I suggest that you instead download the data for chrX from dbSNP's FTP site (this assumes your species of interest is human): ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/chr_rpts/chr_X.txt.gz

Once you have this file, you can use it to filter your dataset:

Collect all rsIDs from the file:

zcat chr_X.txt.gz | awk '/^[0-9]/{print "rs"$1}' > chrX.rsIDs

Filter your dataset.

You haven't shown the format of your dataset, so if it is just a list of rsIDs like this:

$ cat dataset
rs77669712
rs199706910
rs2799940
rs12205639
rs10958205
rs200778574

You can remove any that are on chrX with:

$ grep -vxFf chrX.rsIDs dataset
rs77669712
rs199706910
rs12205639
rs10958205
rs200778574

That will print all rsIDs from the text file dataset that are not also in the text file chrX.rsIDs. Note the -x option for grep which tells it to only consider matches that cover the entire line. This ensures that rs1234 is not removed when you're looking for rs123.


Here a little way to do it. I don't think it is a perfect solution but it should work.
The trick here is to compare patients_id that are wrote with a . or with a -.
So, here is how I do (I modified your example in order to reflect several configurations):

# Reproducible example:
patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
ETV7 = c(1.79,5.38,5,6)
PLAUR = c(6.24,7.05,5,6)
SLC11A1 = c(3.7,3.12,5,6)

df = data.frame(ETV7,PLAUR, SLC11A1)
rownames(df) = patient_id

# OS data
Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
os_data = data.frame(Patient_ids,OS_time)

Now, I create a new column in df that will contains patient_id wrote with a -:

library(rebus.base)
df$Patient_id = gsub(DOT,"-",rownames(df))

Now, I'm looking for common pattern between df$Patient_id and os_data$Patient_ids, find their matching row position and extract these lines:

position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
position = unlist(position)
os_data$Position_matched = position

Now, we are creating a matrix m that we will fill with matching values

m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
colnames(m) = colnames(df)

for(i in 1:length(position))
{
  if(is.na(position[i])){}
  else{
    m[i,] = df[position[i],]
  }
}

And now, we can bind both os_data and m:

New_df = cbind(os_data,m)

And here the final output:

> New_df
   Patient_ids OS_time Position_matched ETV7 PLAUR SLC11A1                   Patient_id
1 TCGA-IG-A4QS    3.93               NA   NA    NA      NA                         <NA>
2 TCGA-JY-A6FB   61.23                1 1.79  6.24    3.70 TCGA-JY-A6FB-01A-11R-A336-31
3 TCGA-JY-A938   35.33                2 5.38  7.05    3.12 TCGA-JY-A938-01A-11R-A37I-31
4 TCGA-JY-A939   22.00               NA   NA    NA      NA                         <NA>
5 TCGA-JY-A93C   23.50               NA   NA    NA      NA                         <NA>
6 TCGA-JY-A93E   25.56               NA   NA    NA      NA                         <NA>
7 TCGA-L5-A4OJ   21.30               NA   NA    NA      NA                         <NA>
8 TCGA-L5-A4OS   59.40                3 5.00  5.00    5.00 TCGA-L5-A4OS-01A-11R-A212-31

Warning: This script should work for some simple case like yours but take extra cautious to the grep output as the short patient_id term of the os_data could eventually generate several hits instead of one.


Demographic model for admixed African Americans

Just to further my answer, the problem with molecular dating within a species (I don't do humans [excepting immunology], but it's the same theory), is it requires an assumption of clonality (bifurcating trees). Mendelian species recombine towards panmixia, which will make any attempt to date via de novo point mutation, appear far younger than they really are, or for panmixia, incalculable (it breaks the central assumption). Dating through admixture requires a lot of extra data and/or assumptions for example the generation time,  putative admixture per generation, and shifts in population dynamics. Maybe you can calculate some of those requirements but  alot will remain assumptions.

I have not read the paper you suggested, because I'm not involved in human genetics [outside immunology], but it is possible the authors may (and I stress may) have simplified an assumption that could be difficult to estimate. It is not clear what the objective of the paper was.

Quantitative analysis is extremely important in modern population genetics/ evolutionary theory because historically it was overlooked and resulted in lots of hypothese which were difficult to test.


Is it possible to check if a patient has the HLA-B27 antigen from his exome stored in a VCF file

Some papers claim that you can do HLA typing from exome sequencing data. I have not validated this myself. I also guess it depends on capture. Nonetheless, I don't think you can type HLA with VCF generated with a standard pipeline, which is not specifically designed for HLA. Current variant callers won't work on HLA genes due to their high diversity. Also in general, there are thousands of HLA-B alleles. In most cases, you can't find the allele by looking at only one SNP.

No, you can't check if the patient has HLA-B27 only from VCF. There might be a hope if you get raw reads and run a proper HLA typing tool designed for such data.


This is what is known as a feature, not a bug. Note that your identical proteins all have accessions starting with WP_. These are special "non-redundant proteins". Many sequences—particularly bacterial sequences—are identical between various different species so having a separate RefSeq entry for each of them would be inefficient. Therefore, RefSeq combines multiple such proteins into a single WP_ record. This is documented here (emphasis mine):

A new type of RefSeq protein record which represents non-redundant protein sequences was introduced in mid-2013. This record type was introduced to address a growing issue with redundancy in the Prokaryotic RefSeq protein dataset that coincided with a significant increase in bacterial genome submissions from individual isolates and closely related bacterial strains. For example, a large number of high-quality bacterial genomes may be submitted during a disease outbreak. The submitted sequences may reflect pathogen evolution during the course of the outbreak but the majority of the encoded proteins from these genomes may be identical to each other. As RefSeq includes these genomes, per community requests, this resulted in increased redundancy. By representing identical proteins using a single non-redundant protein accession number (with the prefix 'WP_'), redundancy in the database is significantly reduced.
[ . . . ]
Because a non-redundant protein sequence may be found in RefSeq genomes from multiple species, the organism information provided on the protein record reflects the lowest-common taxonomic node ranging from the genus species level to super-kingdom. A non-redundant protein record that provides organism information at the level of a genus, family, or even super-kingdom does not mean that the protein is found in all RefSeq genomes below that taxonomic classification. It only indicates that the protein is found in more than one genome of different species for which the genus, family, or super-kingdom classification is the lowest common taxonomic node.

So, your query was NP_418578.1, anaerobic fumarate reductase catalytic and NAD/flavoprotein subunit from E. coli strain K-12, substrain MG1655. The first thing to notice is how specific that is. This is the protein found from one specific substrain of one specific strain of one specific bacterial species. It is reasonable to expect that there will be identical sequences from many, many closely related species. Both from, most probably, all other strains and substrains of E. coli and from other, related bacteria.
Now, the specific sequences you mention are actually slightly different and not 100% identical. Below is a multiple alignment of NP_418578.1 and the 4 WP_ sequences you mentioned. Note that each of the 5 entries is slightly different. Each has one residue that differs from the rest. Look for the : in the identity line, there are 4 : and all others are * (I am only showing the relevant alignment blocks here and have removed those where all 4 sequences were identical):
WP_001401474.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_062863447.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_064226696.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_078165098.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAIAQDH
                    *******************************************************:****

[ . . . ]

WP_001401474.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_062863447.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_064226696.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHIRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_078165098.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ***********************************************:************

[ . . . ]

WP_001401474.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_062863447.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKISQAFWHEWR
WP_064226696.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_078165098.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    **************************************************:*********

WP_001401474.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGIDPVKEPIPVRPTAHYTMG
WP_062863447.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_064226696.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_078165098.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    *****************************************:******************

[ . . . ]

Your sequence (NP_418578.1) is only identical to one WP_* multi-species sequence, WP_001192973:
WP_001192973.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
                    ************************************************************

WP_001192973.1      DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
NP_418578.1         DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
                    ************************************************************

WP_001192973.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ************************************************************

WP_001192973.1      LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
NP_418578.1         LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
                    ************************************************************

WP_001192973.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    ************************************************************

WP_001192973.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    ************************************************************

WP_001192973.1      GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
NP_418578.1         GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
                    ************************************************************

WP_001192973.1      GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
NP_418578.1         GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
                    ************************************************************

WP_001192973.1      AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
NP_418578.1         AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
                    ************************************************************

WP_001192973.1      TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
NP_418578.1         TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
                    ************************************************************

WP_001192973.1      NG
NP_418578.1         NG
                    **

So, in summary, RefSeq will combine multiple identical sequences into a single WP_* multi-species accession. You should therefore expect to find one 100% identical WP_* sequence for your query and multiple, almost identical WP_* entries. And that's precisely what you see here.


How can I convert SCF trace files to ABI files?

Apparently convert_trace does not do a good parameter checking and silently sets ZTR as default if it does not recognise a valid output format.

I don't know of any freely available command-line converters in that specific direction. Even at the height of Sanger sequencing (late 90s, early 2000s), pipelines usually tried to reduce complexity and footprint by converting AB1 to anything else, often SCF (ZTR came too late and never got enough traction).

That being said, I dug out an SCF I had from 1997, loaded it into Geneious, and re-exported it as .ab1. I was a bit astonished to see that this worked. So Geneious (or similar packages from other companies) may be your best bet for a quick, one-off thing.

Depending on your needs, you might want to think about other ways: either expand BioPython to read SCFs, or maybe to expand the Staden io_lib to be able to write .ab1 (James still maintains that package and is extremely helpful).

Good luck.


Why total RNA-seq usually yields low mapping rate?

A likely explanation is that total RNA-Seq contains a high fraction of reads from ribosomal RNAs. Ribosomal RNAs are present in multiple copies across the genome, hence many reads map to multiple genomic locations and get discarded by the aligner. For example, STAR with default parameters considers a read as unmapped if it maps to more than 10 genomic loci (this behavior can be changed with the --outFilterMultimapNmax option).
To confirm if this is actually the case, you could check the number of multi-mapping reads in the aligner's log files.


First of all you are getting an error because you are missing a then in the third line of your code and then a fi to close the if. It should be:

for i in *.vcf
do if grep -q $i 1.txt; then
    cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
fi
done

And I guess the one-liner below achieves what you are trying to do:

while IFS= read -r line; do mv "$line"*.vcf some-directory; done < file-list.txt


With metagenomes, there's several strategies and no one-size-fits-all.
 You can spend a lot of time playing with the data to try to get the "best assembly."  But it depends on what you're trying to do, e.g. whether or not you are mostly interested in the highest abundant species or are interested in the less abundant ones too.

But to propose a strategy that seems reasonable and hasn't been mentioned yet, my initial thoughts are that since the gut microbiome isn't complex, you can likely get decent assemblies from each run separately.  Then you can run a contig clustering algorithm (e.g. MetaBAT) to bin the contigs into genomes.  Next compare the bins from the various assemblies (i.e. cluster the clusters).

It may or may not be worth the trouble to try to combine the related bins' contigs into a consensus (pan)genome sequence.  If you want to, methods that have been used successfully include shredding the contigs into overlapping chunks and assembling with an OLC assembler.  But if you're interested in KO:read-counts (for example), a consensus genome sequence isn't really required.  Cross-mapping the unassembled reads can be useful too (i.e. a sample may have a low-abundance OTU that doesn't produce good contigs, but it can be associated with an OTU that is better assembled in another sample).

If you provide more detail on your goal, I can comment further.


How to best detect the "peaks" in RNA-seq data that are not assigned to any gene?

Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use stringTie or similar programs. Ensure that the BAM file you give to stringTie have all of the cells together.


How to get PSSM of protein by protein sequence with python?

We solve our problem and create gitHub repository (bio-protein).


Is there a point in recalibration of scores for variant calling?

I personally don't think BQSR has a huge impact on variant calling, but you don't really need to guess. If you run GATK BQSR, it outputs a table and charts of exactly how much quality scores are adjusted. The adjustment will vary depending on the position in the read and genomic context (previous and following base). In my experience, the difference is a few points at most, but it's certainly noticeable.

GATK recommends BQSR for both genome and exome data, which is normally much higher than 15x.


Extracting strand specific reads from MinION cDNA-PCR protocol

This has been set up as a bioinformatics protocol on protocols.io.

I've been using LAST to identify the adapter orientation relative to the genome, and then using that information to split BAM files up and recombine them to create two strand-specific files that are displayable in a genome browser.

I have written my own script to process LAST results into a CSV format, which makes it easier to do line-by-line data filtering.

Note: these scripts have been slightly modified from scripts that I have run. Consider them demonstrative: pay attention to the words rather than the script.

Read correction

I prefer starting off my data analysis with a read correction with Canu (ideally v1.7 now that it's out, because that attempts correction of all reads, but here I use canu v1.6). I use minimap as the mapper to speed this up. The genomeSize parameter should be approximately a tenth to a fortieth of the number of bases in your dataset to make sure that no sequences are excluded (bigger is better, as long as Canu doesn't freak out about memory consumption):

~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
  genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
  -p 4T1_BC06 -d 4T1_BC06 -nanopore-raw \
  workspace/pass/barcode06/fastq_runid_*.fastq

This creates a file 4T1_BC06/4T1_BC06.correctedReads.fasta.gz.

Chimeric read filtering

The next step I carry out is a basic read-level QC to exclude chimeric reads. Porechop can be used for this, although that removes adapters by default, which is not particularly useful in this case.

I use LAST to search for adapter sequences within the corrected reads, pass it through my conversion script, and extract out duplicated mappings (i.e. where the same read/adapter pair appears more than once in the mapping results):

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | awk -F',' '{print $1,$2}' | sort | \
  uniq -d | awk '{print $1}' | uniq > reads_with_duplicated_adapters.txt

I use this file to filter out chimeric reads from the corrected dataset using another fastq filtering script I've created:

~/scripts/fastx-fetch.pl -v -i reads_with_duplicated_adapters.txt 4T1_BC06/4T1_BC06.correctedReads.fasta.gz | \
  gzip > 4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Adapter filtering

Along roughly similar lines to the chimeric read filtering, I then look for the strand switch and VNT adapters in the sequences. For a forward-orientated query, I expect the strand switch primer to be in the forward direction, and the VNP primer to be in the reverse direction:

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \
  awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt
lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \
  awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

And then more filtering to split the reads up into forward and reverse [genome-relative] subsets:

~/scripts/fastx-fetch.pl -i fwdQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz
~/scripts/fastx-fetch.pl -i revQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Mapping

Now that the reads have been oriented, the mapping can be done:

~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > fwd_4T1_BC06.CU_vs_mmus.bam
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > rev_4T1_BC06.CU_vs_mmus.bam

Splitting and recombining

The mapped files are then split based on their mapping direction using samtools view flag filtering to exclude or include reverse-mapped reads:

samtools view -F 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > fwd_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > rev_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -F 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > rev_rev_4T1_BC06.CU_vs_mmus.bam

And finally the files are recombined to identify forward and reverse-encoded transcripts via samtools merge. I'm calling the "pos" strand the one that is encoded in the same direction as the genome, and the "neg" strand the reverse-complement direction:

samtools merge pos_4T1_BC06.CU_vs_mmus.bam \
  rev_fwd_4T1_BC06.CU_vs_mmus.bam fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools merge neg_4T1_BC06.CU_vs_mmus.bam \
  fwd_fwd_4T1_BC06.CU_vs_mmus.bam rev_rev_4T1_BC06.CU_vs_mmus.bam

Alternatively, it's possible to set the BAM flags for read 1 and read 2 for the reverse and forward-encoded adapters respectively, which allows everything to then be combined into a single BAM file (and treated in the same way as a strand-specific Illumina BAM file).

If this has worked properly, then mapping human or mouse to the mitochondrial genome should show most expression appearing on the positive strand, with a small scattering of negative-strand expression, a bit like this:

For a bit of a background on how/why this works, see this answer.


Evaluate clusters of individuals by using their sequence data

You can calculate allele frequencies for each cluster you have to further verify if they belong to similar population, however if size of your dataset is rather small this may not work for you.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3813878/, here is an article about differentiating populations based on sequencing data.


Converting aligned fasta to plink ped/bed

An alignment can be the result of two slightly different analyses:
There is multiple sequence alignment (which is what you get from MAFFT) where sequences are aligned so that similar regions are on top of each other. This may require introducing indels (insertions / deletions) if a particular region is absent in some of the sequences.
Then there is alignment to a reference sequence. This is usually how you get variant data and most of the tools deal with this use case. Therefore, in the link you provided the solution includes aligning the data to a reference sequence. Usually the refernce sequence is much larger (e.g. the human genome) and after alignment the result is a file (.bam) that tells you where the query sequences match in the reference.
I did a quick search for ways of converting a multiple sequence alignment (MSA) to VCF. There is a tool called msa2vcf in the Jvarkit collection of utilities that can do this. The example is for the CLUSTAW format but FASTA is accepted as well.
You shoud check that indels were treated correctly because these are most likely to cause trouble.
Then you can simply convert the VCF to PLINK format  using PLINK.


RNA-seq analysis of mixed viral/host reads with salmon

You provide Salmon with a transcriptome fasta... so merging the human transcriptome with the pox virus genome fasta file should work. You don't mention that in your post but the viral genome should exist as a fasta. You can extract the sequences from genbank or maybe you can find our virus here: https://www.ebi.ac.uk/genomes/virus.html

You can find the coding transcriptome at the link below in the cdna folder.
The non-coding transcriptome is in the ncrna. You should probably merge them:

ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/


Intersection of two genomic ranges to keep metadata

Presuming you don't have cases of multiple ovelaps

> library(GenomicRanges)
> gr1 = GRanges(c("chrI", "chrI"), IRanges(start=c(1, 101), width=c(100, 100)), mcols=data.frame(foo=c(0.1, 0.2)))
> gr2 = GRanges(c("chrI"), IRanges(start=c(50), width=c(101)))
> o = findOverlaps(gr1, gr2)
> grl1 = split(gr1[queryHits(o)], 1:length(o)) # You can't mendoapply on a GRanges object
> grl2 = split(gr2[subjectHits(o)], 1:length(o))
> foo = function(x, y) {
+     rv = x
+     start(rv) = max(start(x), start(y))
+     end(rv) = min(end(x), end(y))
+     return(rv)
+ }
> unlist(mendoapply(foo, grl1, y=grl2))

The output is:

GRanges object with 2 ranges and 1 metadata column:
    seqnames     ranges strand | mcols.foo
       <Rle>  <IRanges>  <Rle> | <numeric>
  1     chrI [ 50, 100]      * |       0.1
  2     chrI [101, 150]      * |       0.2
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

It'd be more memory efficient to iterate over o, but this gives you the idea.

In cases of multiple overlaps, one could use intersect(), then findOverlaps() with that and then continue on with the procedure demonstrated above (or iterate over the overlaps).


how is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation?

Note that DIN (DNA Integrity Number) and RIN (RNA Integrity Number) are different scores. I've had trouble in the past finding code or formulas for calculating both of these scores, which is very frustrating considering how frequently they are used in research papers and for NGS QC. The closest that I have been able to get to this are demonstrative graphs for various RIN numbers and DIN numbers, with a description of the variables that are included in the model.

The code that @719016 provided appears to attempt to calculate the RIN, not the DIN.

Just in case it's helpful, here is an extract of that code, only including lines that seem to be directly involved in calculating the actual B_RIN values. This is not runnable, because it depends on particular input files, and it's still pretty unintelligible:

myfiles<-list.files()
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str3<-paste0(substr(targets2[1,2],1,67),"_S")
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
    ##mod chip numb
    targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
    str3<-paste0(substr(targets2[1,2],1,67),"_S")
    samp_num_chip<-length(grep(str3,myfiles))
    ##mod sample numb
    str1<-substr(targets2[1,2],1,67)
    for(j in 1:samp_num_chip)
    {
        iii<-j
        str5<-paste0(str1,"_Sample",iii,".csv")
        samples1[j,i]<-str5
    }
}
samples1<-samples1[grep("2100",samples1)]
##RIN Calculator#######refined peakfinder####################
##read bioanalyzer data into a matrix called dta
##since the total RNA and mRNA assays run differently
##skip more lines for the total RNA assay
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
    x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
    dta[,i] <- x[,2]
}
time<- x[,1]
for(samples in 1:length(samples1))
{
    ##store max fluor in max.peaks matrix
    lad<-samples
    ti1<-17
    ti2<-65
    p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
    peaks <- matrix(NA, nrow=p.row, ncol=1)
    max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
    for (i in 1:length(seq(ti1,ti2,.05)))
    {
        ii<-seq(ti1,ti2,.05)[i]
        x0<- ii
        x0<-round(x0,digits=3)
        y0<- ii
        x<- ii+.05
        x<-round(x,digits=3)
        y<- ii
        max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
    }
    ##human or mouse
    cent<-max(max.peaks)*.6
    ##flatworm
    ##cent<-max(max.peaks)*.05
    ##define the middle of the segments
    max.avg<-seq(ti1+.025,ti2+.025,.05)
    ##find time associated with peaks > 5 fu and fill seconds matrix with them
    ##12.5% of max fu incase chip is fu
    sizes <- which(max.peaks>cent)
    gt35<-which(max.avg[sizes]>35)
    lt55<-which(max.avg[sizes]<55)
    tfr<-which(duplicated(c(gt35,lt55))==TRUE)
    area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],
                max.peaks[sizes][c(gt35,lt55)[tfr]],
                max.peaks[sizes+1][c(gt35,lt55)[tfr]])
    time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],
                 max.avg[sizes][c(gt35,lt55)[tfr]],
                 max.avg[sizes+1][c(gt35,lt55)[tfr]])
    t1<-time2[which(area[,1]<cent),1][1]
    t2<-time2[which(area[,1]<cent),1][2]
    t3<-time2[which(area[,3]<cent),3][1]
    t4<-time2[which(area[,3]<cent),3][2]
    s18a<-ifelse(is.na(t1,-30,round(t1,digits=1)))
    s18b<-ifelse(is.na(t3,-35.05,round(t3,digits=1)))
    s28a<-ifelse(is.na(t2,-30,round(t2,digits=1)))
    s28b<-ifelse(is.na(t4,-30.05,round(t4,digits=1)))
    eta<-which(time==s18a)
    etb<-which(time==s18b)
    tea<-which(time==s28a)
    teb<-which(time==s28b)
    s18<-sum(dta[eta:etb,lad])
    s28<-sum(dta[tea:teb,lad])
    qc.mat[samples]<-  (-1*exp((s28/s18)*-1)+1)*10
}

Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.


GROMACS Coordinates and velocity and force trajectories do not match accuracy

From GROMACS documentation about reproducibility:

The following factors affect the reproducibility of a simulation, and
thus its output:

Precision (mixed / double) with double giving “better” reproducibility.
Number of cores, due to different order in which forces are accumulated. For instance (a+b)+c is not necessarily binary identical
to a+(b+c) in floating-point arithmetic.
Type of processors. Even within the same processor family there can be slight differences.
Optimization level when compiling.
Optimizations at run time: e.g. the FFTW library that is typically used for fast Fourier transforms determines at startup which version
of their algorithms is fastest, and uses that for the remainder of the
calculations. Since the speed estimate is not deterministic, the
results may vary from run to run.
Dynamic linking to different versions of shared libraries (e.g. for FFTs)
Dynamic load balancing, since particles are redistributed to processors based on elapsed wallclock time, which will lead to (a+b)+c
!= a+(b+c) issues as above
Number of PME-only ranks (for parallel PME simulations)
MPI reductions typically do not guarantee the order of the operations, and so the absence of associativity for floating-point
arithmetic means the result of a reduction depends on the order
actually chosen
On GPUs, the reduction of e.g. non-bonded forces has a non-deterministic summation order, so any fast implementation is
non-reprodudible by design.

I would add that in your case of hard coding everything by hand you may have to look a bit into your barostat and thermostat (although their influence should be found on speed of the atoms). Maybe a sample code with coordinates and speed of 4 atoms would help to give a more precise answer.
You can also try having a look at gromacs source code. I believe OpenMM provides a bit more control and clarity on each step, if for some unfathomable reason you are not doing this only for curiosity sake but because you actually have to make a code that reproduces these results.


Missing genes and normalisation of RSEM output using EBSeq

Yes, that blog post does represent just one guy's opinion (hi!) and it does date all the way back to 2014, which is, like, decades in genomics years. :-) By the way, there is quite a bit of literature discussing the improvements that expected read counts derived from an Expectation Maximization algorithm provide over raw read counts. I'd suggest reading the RSEM papers for a start[1][2].

But your main question is about the mechanics of running RSEM and EBSeq. First, RSEM was written explicitly to be compatible with EBSeq, so I'd be very surprised if it does not work correctly out-of-the-box. Second, EBSeq's MedianNorm function worked very well in my experience for normalizing the library counts. Along those lines, the blog you mentioned above has another post that you may find useful.

But all joking aside, these tools are indeed dated. Alignment-free RNA-Seq tools provide orders-of-magnitude improvements in runtime over the older alignment-based alternatives, with comparable accuracy. Sailfish was the first in a growing list of tools that now includes Salmon and Kallisto. When starting a new analysis from scratch (i.e. if you ever get the original FASTQ files), there's really no good reason not to estimate expression using these much faster tools, followed by a differential expression analysis with DESeq2, edgeR, or sleuth.

1Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN (2010) RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics, 26(4):493–500, doi:10.1093/bioinformatics/btp692.

2Li B, Dewey C (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12:323, doi:10.1186/1471-2105-12-323.


If you want to find the intersection between ALL your bedfiles, you can try multiIntersectBed (available since bedtools 2.14.3).

It should work like this :

multiIntersectBed -i SRR2920506.filtered.bed SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

To sort your bedfiles by chromosome then by start :

sort -k 1,1 -k2,2n myfile.bed


Samtools Index: Chromosome Blocks not Continuous

samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

These error messages indicate that the reads are not sorted by coordinate — in particular, that the reads mapped to ScdB1pO_646;HRSCAF=880 are not all together, and that the reads at positions 11159020 < 11158717 are not sorted by position on their chromosome.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.


cn.MOPS fails with 'missing value where TRUE/FALSE needed'

I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome chrUn_gl000228. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the exomecn.mops(X) was failing.

So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on chrUn_gl000228 and confirmed that the program worked as expected and didn't throw any errors.


Getting protein FASTA sequence based on keyword with python

Found what I was looking for. In:

searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)

I've added:

searchterm = "(terminase large subunit AND viruses[Organism]) AND Caudovirales AND refseq[Filter]"
searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=6000)

which norrowed down my searches to the desired viruses. Granted it's not filtered by host, but by a taxonomy group, but it is enough for my work.

Thank you @Llopis for additional help


tSNE often offers better visual representation (separation) on such complicated data than PCA. As Micheal pointed out, computing a tSNE embedding over 20.000 gene dimensions is computationally unfeasible, so a number of PCs are normally calculated and these are used as input for calculating the tSNE. They are used in tandem.

As for global vs. local, we are much more interested to see similarity of cells to a limited number of neighbours indicating a celltype and grouping these close together. This is more important than the distance between such celltypes. (assuming the separation in your tSNE is driven by a biologically meaningful factor such as celltype    and not some confounder)

Edit:

Since I just made these images for myself anyway I might as well post them. The first 2 PCs show some separation by celltype, tSNE computed over 100 PCs gets very nice separation.


Filtering imputed GWAS SNPs based on a MAF difference of 10%

I cannot think of any principled rationale for choosing this filtering strategy.

However, I am going to take a guess that this filtering strategy is supposed to filter out SNPs for which imputation did not work well?  In that case the appropriate statistic to filter on is the INFO score as described here.

You might consider picking a higher threshold than 0.4 if you feel that you aren't filtering out enough variants.


Trim Seurat object to contain expression info only for selected genes

A different approach if you are using Seurat3, is DietSeurat(). It allows you to diet the object by removing the components that you don't need. For example you can keep the normalised/scaled matrix and remove the raw counts. This approach could reduce space and memory usage, while keeping all your genes in place.


Visualizing the difference among three groups

Here a way to do it is to start first by creating a dataset containing your three different group, based on the few lines you display:

dat1$Condition = "dat1"
dat2$Condition = "dat2"
dat3$Condition = "dat3"
colnames(dat1)[1] = "Patient"
DF <- rbind(dat1,dat2,dat3)

             Patient   DEL   INS    SNP  total Condition
1: LP6008337-DNA_H06   927   773  40756  42456      dat1
2: LP6008334-DNA_D02  1049   799  31009  32857      dat1
3: LP6008031-DNA_E01 13552  3374  62105  79031      dat2
4: LP6005500-DNA_G01   539   500  43451  44490      dat2
5: LP6005935-DNA_F03 39168 16739  58095 114002      dat3
6: LP6008269-DNA_D08   849   910 103501 105271      dat3

Then, you can reshape DF in order to put all values DEL, INS, SNP and Total into a single column and create a categorical column. Doing this will be helpful for later producing the facets of your plot. To do that, you can use pivot_longer from tidyr package:

library(tidyr)
DF <- DF %>% pivot_longer(.,cols = c(DEL,INS,SNP, total), names_to = "var", values_to = "val")

# A tibble: 24 x 4
   Patient           Condition var     val
   <chr>             <chr>     <chr> <int>
 1 LP6008337-DNA_H06 dat1      DEL     927
 2 LP6008337-DNA_H06 dat1      INS     773
 3 LP6008337-DNA_H06 dat1      SNP   40756
 4 LP6008337-DNA_H06 dat1      total 42456
 5 LP6008334-DNA_D02 dat1      DEL    1049
 6 LP6008334-DNA_D02 dat1      INS     799
 7 LP6008334-DNA_D02 dat1      SNP   31009
 8 LP6008334-DNA_D02 dat1      total 32857
 9 LP6008031-DNA_E01 dat2      DEL   13552
10 LP6008031-DNA_E01 dat2      INS    3374
# … with 14 more rows

Now, you can plot your data and use the function stat_compare_means from ggpubr package to get the significance display on each tab. Based on the size of your group, you may want to adjust for the type of statistical test to be used (take a look at ?stat_compare_means to find the statistical test adequate to your dataset). You can pass various conditions to test as a list by using comparisons argument.

For facetting your graph, you can use facet_grid (orfacet_wrap, try the one you preferred). Altogether, you can get the plot by writing:

library(ggplot2)
library(ggpubr)
ggplot(DF, aes(x = Condition, y = val, fill = Condition))+
  geom_boxplot()+
  stat_compare_means(comparisons = list(c("dat1","dat2"), c("dat1","dat3"), c("dat2","dat3")))+
  stat_compare_means()+
  facet_grid(var~., scales = "free")

Hope it helps you to figure it out how to represent your data.


Since every human has a different DNA (different combinations of C, G, A, T) what does it mean to have the genome done?

All humans have some differences in their DNA, but there's far more that is shared. On average the difference between humans is only about one thousandth of their full DNA, which means we're about 99.9% the same. These differences aren't distributed fully randomly, but are often because of specific gene alternatives. (Random mutations do occur, but they are also often fatal, so the random mutations we see of living adults are much more restrictive than all the random mutations that occur.)

To identify the human genome is to study many people's DNA and to label the parts that are shared between everyone, the parts with two or three variations, and the parts with even more variations. Even though every person has different DNA, we can still say they fit the pattern, just like every T-shirt may be unique but they all fit the T-shirt pattern and not the trousers pattern.


NO_COOR reads not in a single block at the end 0 -1

This is essentially a duplicate of Samtools Index: Chromosome Blocks not Continuous but with different error messages, and a duplicate of Error given while trying to index a BAM file with Samtools Index - NO COOR?.
samtools sort -n ${file} -o ${file}_sorted.bam

These error messages also indicate that the reads are not sorted by coordinate — in particular, that the unplaced unmapped reads with RNAME * are not all together at the end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.


Why are TPMs per 10k or 100k in many scRNA-seq studies?

Collecting all the question comments as a community answer (feel free to edit this to make it more readable):
There tend to be far fewer counts in single cell experiments. By reducing the scaling factor you are avoiding having massive expression values for your samples.
It isn't uncommon to scale single cell experiments by the median cell library size. That way small values typically won't be too small and large values won't be too large. People also choose 10,000 because that's roughly around the library size for a single cell.
I haven't seen an explanation yet. Here is an example paper from last week (reputable journal and lab): doi.org/10.1016/j.celrep.2018.11.056 . "... We also excluded outlier cells from the MPPs group which were likely caused by cell contamination. To correct for batch issues, the two experimental batches were centered to have the same mean log(TP100K+1) (hereafter referred to as log(TPM), where TPM stands for ‘transcripts per million’) per gene. PCA was done in R with scaled..."


Here's how I'd approach this:

for READ1 in /Documents/Data/Input/*R1*
do
    READ2=${READ1/R1/R2}
    mixcr analyze amplicon --s hs \
        --starting-material dna \
        --5-end v-primers \
        --3-end j-primers \
        --adapters no-adapters \
        --receptor-type trb \
        --only-productive "$READ1" "$READ2"
done

Basically, loop over all the R1 FASTQ files in a directory, and calculate the R2 names as needed.

This is using BASH string manipulation, and assumes that substituting R2 for R1 will correctly calculate the R2 filename from R1. It also assumes that all the FASTQ files are in a single directory.


Find gene at position from gff or gbk file

Via BEDOPS:

$ gff2bed < annotations.gff > annotations.bed
$ vcf2bed < snps.vcf > snps.bed
$ bedmap --echo --echo-map-id-uniq snps.bed annotations.bed > answer.bed

This can be reduced to a one-liner if you're using bash:

$ bedmap --echo --echo-map-id-uniq <(vcf2bed < snps.vcf) <(gff2bed < annotations.gff) > answer.bed


Is there a good practice (or an easy one for a SBOL newbie) to set up a collection of biological parts and devices?

For SBOL 2, you can readily import data via the online converter tool at https://converter.sbolstandard.org/
The converter basically just wraps the conversion function found in the SBOL libraries, so if you're coding, you can do the conversion directly from the libraries.
For SBOL 3, there isn't yet an online tool, but James McLaughlin has a converter in his JavaScript library: https://github.com/sboltools/sbolgraph
I believe he is working to bring it online shortly, so hopefully this answer will soon need to be updated or replaced by an answer from him.
An anecdote on the difficulty level: I recently started using an SBOL3-->SBOL2-->GenBank conversion using pySBOL3 and pySBOL3. It took a couple of hours to setup the pipeline, and isn't perfect, but it's pretty doable (and ought to be generalized for others to use as well at some point). The other way would be easier, since GenBank contents are pretty simple.


How to calculate statistical significance of sequence motifs

As far as I remember the exact probability computation is an open problem. The reason is that potential motifs can overlap, which makes probability computations for an arbitrary string non-trivial, and depends on the motif.

For example if he have a binary string of four digits, the probability of "01" will be 11/16, while the probability of "11" will be 8/16=1/2.

The simplest approximation is to assume that probability of the motif on every position of the sequence is equal and independent. In such case (this comes from this course):

(Probability of text of length N having at least t occurrences of a k-mer Pattern, A is the number of letters in the alphabet, n = N – t * k)

This approximation is very rough, there are better here.


How to perform functional analysis on a gene list in R?

N.B., I'm avoiding discussion of "best", since that's more or less impossible to answer.

Your question can actually be divided into two:

What's a good tool for pathways analysis (ideally in R)?
What are good sources of pathway information?

For (1), there are a few different possibilities, but I prefer either roast or camera from limma or goseq (a stand-alone bioconductor package). Of those, I expect what roast or goseq are doing are the most similar to what panther is doing in the sense that they're not competitive. However, since genes do correlate with each other, I prefer competitive tests like that offered by camera.

For (2), your options are basically KEGG or reactome. KEGG is problematic, since you need a license for anything remotely recent (there are some R packages to get around this, but I haven't a clue how kosher they are legally). Given that, reactome might be your best bet.

Having written all of that, the best database out there is what IPA uses. This is a commercial product (I have no affiliation), but you can get a demo license that will work for a couple weeks at least. If you're doing analyses relatively often then it makes sense to spend a bit of cash on a license (maybe share with neighboring groups).

Regardless of what you use, you might also want to use pathview, which can produce easily interpretable plots related to pathway enrichment.


Yes, just add up the coverage at each location; that's integration in a nutshell.

Or, alternatively / equivalently, count the total number of reads within the region.


How to apply function on bedfile based on a reference bedfile

You could munge your first file into a BED5 file:

$ awk -vOFS="\t" '{ print $1, $2, $3, $6, $5 }' /tmp/f1.bed | sort-bed - > /tmp/f1.bed5

Sort the second file, if necessary:

$ sort-bed /tmp/f2.unknownSort.bed5 > /tmp/f2.bed5

Take the union of the two sorted files:

$ bedops -u /tmp/f1.bed5 /tmp/f2.bed5 > /tmp/union.bed

Here's what the union looks like:

chr1    1       5       .       0
chr1    1       5       Sample2 1
chr1    1       10      Sample1 1
chr1    5       20      .       1
chr1    20      180     Sample1 1
chr1    50      120     .       -1
chr1    70      100     Sample3 1

We have one name scheme for reference element IDs (.) and a second name scheme for sample IDs (Sample*).

We will exploit this difference in name schemes when we apply operations in a certain order. For this example, we will subtract a reference element score from a sample element score.

We can partition the union, mapping the disjoint elements back to the unioned set:

$ bedops --partition /tmp/union.bed | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed
chr1    1   5   .;Sample2;Sample1   0.000000;1.000000;1.000000
chr1    5   10  Sample1;.   1.000000;1.000000
chr1    10  20  .   1.000000
chr1    20  50  Sample1 1.000000
chr1    50  70  Sample1;.   1.000000;-1.000000
chr1    70  100 Sample1;.;Sample3   1.000000;-1.000000;1.000000
chr1    100 120 Sample1;.   1.000000;-1.000000
chr1    120 180 Sample1 1.000000

Each of the fourth and fifth columns are semi-colon-delimited strings. Each element is ordered by the order in which it is discovered in a bedmap operation.

Both columns have the same ordering or same one-to-one correspondence of ID to signal score. This is very useful, because we can turn these into a pair of ordered arrays and apply some logic, depending on the answer we want back.

For example, here is some awk logic to parse those two strings and subtract a reference element score from a sample element score:

$ bedops --partition /tmp/union.bed \
    | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed \
    | awk -vFS="\t" -vOFS="\t" \
        '{ \
             i=split($4,a,";"); \
             j=split($5,b,";"); \
             if ((i==1) && (a[i]==".")) { pass; } \
             else { \
                 for (k=1; k<=i; k++) { \
                     if (a[k]==".") { refIdx=k; } \
                 } \
                 for (k=1; k<=i; k++) { \
                     if (k!=refIdx) { print $0, ".", b[k]-b[refIdx], a[k]; } \
                 } \
             } \
         }' \
    | cut -f1-3,6-8 \
    | sort-bed - \
    > answer.bed

The file answer.bed will look like this:

chr1    1   5   .   1   Sample1
chr1    1   5   .   1   Sample2
chr1    5   10  .   0   Sample1
chr1    20  50  .   1   Sample1
chr1    50  70  .   2   Sample1
chr1    70  100 .   2   Sample1
chr1    70  100 .   2   Sample3
chr1    100 120 .   2   Sample1
chr1    120 180 .   1   Sample1

This partitioning of intervals is different from the one in your question. Perhaps this is due to how you partitioned your set. In any case, perhaps this would give you an idea of how to apply custom operations over a disjoint partitioning.

The main idea is that this part of the awk script finds the score associated with the reference file ("File2"):

for (k=1; k<=i; k++) {
    if (a[k]==".") { refIdx=k; }
}

Any index other than refIdx corresponds to a Sample* element. That means we can do operations between Sample* scores and the specific reference score.

We also leave out mappings where a reference score does not associate with any Sample* element:

if ((i==1) && (a[i]==".")) { pass; }

Thus, all the operations are always between a pairing of a sample and the mapped reference element.

The cut and sort-bed statements at the end just clean up the result so that you get something that looks like what you are expecting at the end, which is also sorted correctly and ready to use for downstream set operations.


From Wikipedia:
In the mathematical field of dynamical systems, an attractor is a set of numerical values toward which a system tends to evolve, for a wide variety of starting conditions of the system.

In biology you can think of an attractor as a state towards the system is moving to. A good example is the process of differentiation where stem cells compromise to more specialized cells, that could be thought as the attractors of the process.


How to compare the contents of a column of the same data table

Loop through 3 columns by row using apply, then convert toupper case (as we want a same as A), get unique, and get lengths, if length is 1 then all 3 values are equal:

df1[ lengths(
  apply(df1[, 1:3], 1, function(i) unique(toupper(i)))
  ) == 1, ]

This would scale better if we want to compare 4, 10, n(?) columns.


If you only have a small number of genes, OMIM is a very easy way to do this.
For example:

If you have too many to do manually, I suspect there is an easy way to query lists, although I'm not sure exactly the way myself. Perhaps someone can add a new answer if so.


For discordant mapped reads, a position of 0 for the pair of a mapped read usually means that the pair couldn't be mapped anywhere. This would be expected for fusion events because some reads could originate from across the fusion site (i.e. not being fully contained in either "unfused" gene), leading to a sequence that doesn't match anything in the reference.


Purpose of ### (3 consecutive pound signs / hashtags / octothorps) in GFF3

The three # are used for splitting group of features that belong together, e.g. a transcript and it's exons. Sometimes you see with a blank line instead of the #.

To grep just the header, which has 1 or 2 # from the line start on, you can use extended regular expression:

$ grep -E "^#{1,2}[^#]" Caenorhabditis_elegans.WBcel235.95.gff3

Which means: "grep all lines that started with 1 or 2 # and are not followed by #"

If you like to create a new gff3 to get rid of this separator lines you can do this:

$ grep -v "^###" Caenorhabditis_elegans.WBcel235.95.gff3

The -v prints all line that doesn't match the pattern. (I guess this is what @Michael mean by "negator")


Your understanding of FASTA format is about right. The type of basic problem you're eluding to we term "sequence alignment"- edit distance might be okay for teaching but in practise we use other algorithms, e.g. you might be interested in the Needleman–Wunsch or Smith–Waterman algorithms. Richard Durbin et al. wrote a great book that covers these any much more https://www.amazon.com/Biological-Sequence-Analysis-Probabilistic-Proteins/dp/0521629713
Tracking mutations etc. requires more than just alignment though, see "phylogenetics" (i.e. building genetic trees) and "variant calling". Also checkout what the nextstrain team is doing https://nextstrain.org/ncov/global
In general, looking for practical applications for your algorithms is great, but be very careful before drawing any real world conclusions about the coronavirus outbreak from such analyses


Can I use DRIMseq to model interaction effect of variables?

One way to see the interaction effect would be including it in the design matrix as explained in the limma tutorial. Something on the line of:

design_full <- model.matrix(~ treatment*strain, data = samples(d))

Or you could build it yourself in a matrix like

 (Intercept) strainA_treatmentA strainA_treatmentB  strainB_treatmentA strainB_treatmentB
       1                      0                  1                  0                   0
       1                      1                  0                  0                   0

About the error, it seems like the data is not in the right format, you have either NA or Inf somewhere in the matrix. Check it and rethink the processing you made.


The help page for removeBatchEffect() explicitly states that it should not be used in conjunction with differential expression testing, so don't do that. You already have the batch effect included in your model, so you're already correcting for it when testing for differential expression. If your MDS plots suggest that samples are variably affected by the batch effect then look into the combat() function from the sva package.


Are the conclusions in "The proximal origin of SARS-CoV-2" legit?

In summary, the authors are saying the complete opposite of "human intervention".

  While the analyses above suggest that SARS-CoV-2 may bind human ACE2
  with high affinity, computational analyses predict that the
  interaction is not ideal and that the RBD sequence is different from
  those shown in SARS-CoV to be optimal for receptor binding.

The interaction with ACE2 is taking place on the Spike protein (S) in both COVID-19 and SARS-Cov. The S protein is very different between COVID-19 and SARS both in terms of sequence and COVID-19 includes an additional furin cleavage site, meaning structurally it is quite distinct. The "computational analysis predicts" is difficult because there is no crystal structure for COVID-19 S protein therefore modelling ligand binding is essentially guess work, particularly as it comprises an additional (furin) cleaved domain. What the authors appear to conclude is that they cannot map ACE2-Spike binding of COVID-19 using an established SARS S protein structure. I assume they are homology modelling SARS-CoV onto COVID-19 S protein sequence data.  In my personal opinion there is no way to homology model an additional furin protease cleavage site (furin is a host protease) so it would be difficult for "good fit" to occur, but it is just my opinion and the authors may have established methods to overcome this notable limitation, e.g. similar phenomena occur in influenza.

  Thus, the high-affinity binding of the SARS-CoV-2 spike protein to
  human ACE2 is most likely the result of natural selection..

Their conclusion of "natural selection" is not strictly accurate and should be "Darwinian positive selection", but its a small issue. What they are saying is Darwinian adaptation has occurred as opposed to "purifying selection, i.e. conservation. Thus they are definitely excluding human intervention because that is not "natural selection"

  Furthermore, if genetic manipulation had been performed, one of the
  several reverse-genetic systems available for betacoronaviruses would
  probably have been used

They appear to be saying that if a synthetic virus was constructed it would be created using an established reverse genetics system, and all those produced to date look nothing like COVID-19.

The phylogenetics argument is always a stronger argument than either of the above (although my opinion isn't entirely objective).  Phylogeneticists have had a long history in opposing conspiracy theories, e.g. polio vaccination, HIV and the list goes on, so it is much more natural territory. An example of such an argument is here which @terdon helpfully pointed out.

To address the questions of @Hans

Question 1

Purifying selection = any evolutionary change is deleterious, so the virus is less able to transmit between humans. The vast majority of mutations are deleterious in classical Darwinian thinking*
Adaptation = the amino acid change occurs to ensure the virus is better able to transmit between humans. In phylogenetics this is called 'positive selection' which is detected at a nucleotide level, this is however a very stringest test.

*, There is something called nearly neutral theory, but thats just complicated

Question 2

What we are talking about is a virus infectious clone system. For a 30kb virus this would not be trival at all, but I agree for 10kb viruses it much easier. However, the phylogenetics are opposed to this, because SARS-CoV-2 is 5% divergent from a virus (RaTG13) isolated a long time before the current epidemic. You can't just engineer a virus that is 5% divergent from a bat virus across the entire genome and never been seen before. It would be a work of genius the likes we have never seen using bioinformatics as yet unheard of, because what would happen is you would encounter endless deleterious mutations. Even if you got over all those, how do you know what fitness traits you are aiming for?

Serial passage
This level of change is doable by serial passage of the virus, but this leads to attenuation of virus (less dangerous) not increased fitness. The yellow fever virus vaccine was created in this way. Even if a virus had been serially passaged we could tell, because we know the mutational patterns of other serially passaged viruses.


In cases when you cannot find answer through the web interface, you can download all the data and search it locally.

Relying on the PDB annotations of SS bonds, you could search for PDB files that have exactly one SSBOND record:

$ cd $PDB_DIR/structures/divided/pdb/
$ zgrep -c ^SSBOND aa/pdb*.ent.gz | grep :1$
aa/pdb2aaq.ent.gz:1
aa/pdb3aad.ent.gz:1
aa/pdb5aap.ent.gz:1
aa/pdb6aa4.ent.gz:1
aa/pdb6aa5.ent.gz:1

All PDB files (gzipped) take about 30GB of disk space.

The caveat is that some entries are available only in the mmCIF format. Using mmCIF would be slighty more complex: you would need to count _struct_conn.conn_type_id being disulf. And you would need to download 10GB more.


How to draw a 3d structure using 3d coordinates from pdb file?

JSmol sounds like what you want - it'll let you display 3D structures from PDB files onto webpages using JavaScript. If you just want to view the structure rather than display them on webpages, there's many programs that exist already like Swiss-PdbViewer and Pymol. There's even web-interfaces that let you upload a PDB file and they display it for you like Jena3D Viewer.


Compute the significance of the overlap between 2 or more gene sets

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

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

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

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

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

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

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


Existing tool for converting gff3 to genbank (gbk)

I believe seqret is already the simplest approach. If you install the emboss command line utility, genome (file) sizes will not be a problem. You can even install it with conda, see https://anaconda.org/bioconda/emboss.

After that, you can use it as the example:

seqret -sequence {genome file} -feature -fformat gff -fopenfile {gff file} -osformat genbank -osname_outseq {output prefix} -ofdirectory_outseq gbk_file -auto

Hope it helps


Getting exact mismatch data from multiple sequence alignment

I have partly answered this for know! For anyone interested i use kalign to align all 10 strains to a fasta format, then I use snp-sites to output a vcf file that can show each mismatch in accordance to the first strain in the fasta. New problem is that the alignment "stretches" the length of this first strain, because some of the other strains are a bit longer. This means that the position of variants are with reference to the longest strain of the 10 and not the strain I want. It can work for now, but any solution to this, would be helpful.


using snakemake shadow rules to store temp files on local nodes

What is happening is your Admin is not allowing you access to /scratch ... which is undertandable. Essentially you are attempting to access the headnode /scratch. If every user on the cluster had access, it would be a bit chaotic.
I am surprised by the error, not that I would use Snakemake I'd code it, because you would expect a data pipeline tool as established as Snakemake to have a solution to working on a cluster. What you are describing is a standard cluster set-up between headnode and 'sub-node' (the technical term is slavenode, but I hate that term). The headnode is mirrored across the cluster, hence the error. Once you submit via qsub the mirror will 'display' the headnode layout on the subnode.

You could ask your admin for permission to /scratch, this is root privilege. It is a bit unlikely they will say 'yes' but its worth a try.
Otherwise what about a virtual environment? When Perl suffers this problem I simply use 'perlbrew', anaconda has very nice virtual environments that are easy to setup
@Chris_Rand points out piping to a different output.
I bet there is a solution in '.profile' and thats how I solved comparable problems in the past.
Of course Snakemake do a Q&A interaction on Overflow SE (but please come back here!) ... they must know the work around.

Frankly I am astonished because the whole point of Snakemake, is the ease of constructing a data pipeline and data pipelines are processed on clusters, thus this error seems to up-end it, but who am I?


hisat2 --rna-strandness option and downstream htseq-count analysis

If you reran the command with the correct settings, just leave it at that.  (It is not at all clear to me that strandedness rf is correct)

If you want people to tell you if you ran the commands right, you need to put down what commands you used.


What algorithms exist for reverse protein folding?

The researchers behind FoldIt have the occasional design challenges that use human game players to help design proteins to fit a structure.

However, even with a known protein sequence there's still a lot of labwork involved in trying to find a good working DNA sequence. Just generating the DNA sequence is a difficult task due to redundancy in the genetic code (i.e. 64 codons for 20 amino acids + stop signals). When that is combined with cell-specific tRNA frequencies that can influence the timing of amino acid incorporation (and therefore how a protein folds as it is created), creating a sequence that reliably fits a shape ends up being quite difficult.


How to select the most representative pathways from a gene enrichment analysis?

This sounds like something that might be amenable to a clustered heatmap plot, or a correlation matrix plot, or something similar. Have you looked at a correlation matrix of the dice coefficient matrix (or maybe just a heatmap plot of that matrix without the correlation matrix)?

The corrplot package looks like it might be useful, in particular the hclust / drawing rectangles presentation.

I can't vouch for this package though; it's just something I found by a search for "R plot correlation matrix".


Interpreting contig alignments to a reference genome

some additional information to your request.
https://ftp.ncbi.nlm.nih.gov/genomes/Viruses/FamilyPhylogeneticTree/Graph/SRR12707453.MN813694.1.gfa
it contain assembled by reference graph
using ncbi tool saute(skesa) which will be public later.
install bandage app and load this graph
https://rrwick.github.io/Bandage/
picture looks intriguing,I would say , run contain many errors located in one position.
Reason? I don’t really know.


Genetic Genealogy related to SNP variants and INDEL variants

InDels are just as informative, but are less reliable to call. Errors should be averaged out in this kind of analysis though (so long as they don't predominate).

Both SNPs and InDels are valid data points to consider when comparing triads, but identifying IBD segments will only ever be probabilistic.


Evaluating amateur bioinformatics results

You need to understand the biological significance of the 33-mer locus.

I'd blast your 33-mer on NCBI and see what the biological significane is of the locus. If might be a transposon/ promoter region.
Look at the translation of the protein and if its part of a protein the location within the protein structure.
Blast the entire locus including 5' and 3' regions of the 33-mer target
Carefully assess the results for chance, e.g. conserved protein and the synonymous mutations just happened to drop out (Poison distribution).

Once you gain a greater understanding, it is likely to have been identified before and what you might contribute is a different method to obtain a known result.

In summary, you need the next level part of the course about translating a genetic correlation into a biological result/ mechanism. This isn't research grade work at present in my opinion.