text:
filter secondary alignments using pysam
I am new to python and trying to learn. The below is an attempt to filter out secondary reads in a bam using pysam. I ran mark duplicates and got an error on several bam files due to secondary alignments. I added comments to each line to illustrate my thought process.  Thank you :).
The specific mark duplicates error:
Value was put into PairInfoMap more than once

Here is the script:
#! /usr/bin/python  ## call python script
import sys ## import python system functions
import pysam  ## import module

bam = pysam.AlignmentFile(*.bam, "rb")  ## open bam and read
click.echo("Reading BAM file")  ## output message
for read in bam:  ## start loop and iterate over each bam
    if read.is_secondary=true  ## not the primary alignment
    if read.has_tag('XS')  ## look for Secondary Alignment Score tag
        XS = read.get_tag('XS')  ## read and store AS value
    if read.has_tag('AS')  ## look for Alignment Score tag
        AS = read.get_tag('AS')  ## read and store value os XS tag
    if AS < XS  ## Alignment score less then
        read.write(read)  ## only print primary alignments
bam.close()  ## end processing

Is there an issue with the script?  If so what is it?  Are you looking for feedback on the script?

Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the read.has_tag('XS') continue block.
I think that if you filter out is_secondary and is_supplementary you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly.
for read in bam.fetch()
    if read.is_secondary=true
        continue
    if read.has_tag('XS') <- This block stores XS and then moves on to the next read
        XS = read.get_tag('XS')
        continue
    if read.has_tag('AS')
        AS = read.get_tag('AS')
    if AS < XS  <- This block compares the AS of the current read to the XS of the previous read that had one.
        read.write(read)
bam.close()

I have not tried to executec it as this is my first python script i am looking for suggestions. Thank you :).

I would suggest running the script to determine if it is doing what you want.  What I pointed out above is def a flaw in the logic.

I ran the script (modified it in the post) and there are many syntax errors... being new to python I am not sure what is best... what do you suggest? Thanks :).

I would start will researching how to write to a file.  read.write(read) is not correct syntax.  You will need to open a file and write to that.  Go error by error until you have a working script then determine if the out is correct.


text:
Interpretation of hierarchical clustering with bootstrapping
I have data that includes 'cases' and 'controls' and have carried out hierarchical clustering.
I used the pvclust package in R to bootstrap the results and significant branches are highlighted with red rectangles (based on au>0.95):

What is clear is that no clustering occurs that separates 'cases' and 'controls', and this is in fact what we expected and want to show. We want to show that the measured variable does not distinguish between cases and controls.

List itemApart from saying visually no clear clusters emerge that distinguishes between cases and controls are there any objective measures that can be used to say no significant clustering occurs between two groups?
One observation I have here is that the AU values and the BP values are very different, even though both p-values should be interpreted in a similar fashion, am I missing something?
Perhaps pvclust (bootstrapping) is not the right option here, is there a better way of showing quantifiably that no significant clustering occurs between two groups? Perhaps some kind of supervised clustering (I am not sure what this even means in this context)?

I am unfamiliar with this package but some red flags: 1) you use correlation as a distance measure, it makes more sense as a similarity measure, 2) I don't see a cutoff in the tree to determine what your actual clustering is. This might give a view where clusters do emerge with distinct case/control patterns. edit: I see the red brackets on the bottom now, are these your final cluster labels?

Can you explain more about your data please? Is this gene expression? Or something else?

@Pallie, I have cut off a section of the picture which has patient ID's in it so there is a blank white space which obscures the red boxes sorry. If you look at the red boxes on the dendrogram they highlight branches with AU values >95 which I interpret as meaning these are the statistically significant clusters.
Benn, this I not gene expression, it is a measure of antibodies to about 100 bacterial strains in cases vs controls.
Any advice on a better measure to use perhaps to show the measure does not cluster patients and controls? Thank you both

I'm not sure doing clustering and then trying to prove that the clustering is significant is the right approach, it is very heavily dependent on the clustering algorithm you used. Can't you just do an ANOVA between case/control to prove there is a/no significant difference?

ANOVA requires normalised data and its unlikely to be Gaussian (it sounds a bit odd saying "unlikely to be normal).

The take-home lesson from this question is always ask what data it is. I'd assumed it was some sort of gene expression, antibody data is radically different and explains the results better. The problem is that IgG will be capable of binding to some degree with any strain in the analysis.

Thanks, yes lesson learned, sorry for the confusion, allot to learn.

Good question. You have performed the first step via unsupervised learning and AU bootstraps give great results, orthodox bootstraps give no significant (or very few) clusters. The contrast between these approaches is unusual because they should be somewhere close to each other. If you can resolve that, your analysis is fine. So at present IMO the unsupervised learning needs further investigation.
... deleted
... antibodies ... oh I see. Trees, like you've done, don't work here because of the cross-reactivity creates network phylogenies. This is the antonym of bifurcating trees. If someone has a robust bifurcating phylogeny with antibody data and gets good bootstraps .. utterly amazing. However, if you subject your data to a network analysis you'll find connections all over the place - in practical terms this really messes up bootstraps. Thats the likely reason for you results. There is no point doing a complex bifurcating tree analysis because the antibody cross-reactivity will mess it up. k-means has been implemented as unsupervised learning in these instances and would provide the starting point in using clustering to define a supervised learning problem.

could you briefly explain how K-means clustering would work here? would I simple choose k=2 (cases vs controls) and show that cases and controls do not cluster into 2 groups? Could I still bootstrap the result to give significance?

No not in this instance ... with "traditional" k-means is you have to set the groups a priori 2,3,4,5 etc..... and the trick is defining which a priori grouping provides the "best fit" (I forget the criterion). The tree will provide some indication of the max group size. K-means will then assign all you data into the designated number of clusters as best it can. K-means is respectable - albeit it doesn't sound cool. I think its because it minimises false positives and robust. The study I was involved with (malaria) used k-means as the singular proof :-O

In this study, the network phylogeny defined the groupings (it couldn't assign significance to it) and these groupings were explored via K-means. Trees (network trees) were there just because I like them, all the tree was saying that was statically significant was "these things cross-react" (hardly revolutionary). Otherwise it is a bit hard to describe k-means within a comments section, but let me know if you get the gist.

I am doing supervised learning on quite similar data to yours at present but for viruses. However, this is a bit too close my current work to be providing detailed technical advice.

I've got it.
The question you are asking is homogeneity between you target group and your general population of patients. What you want to show is "complete" homogeneity between the two groups regardless of any sub-sampling, ie. if there is groups within groups that might break your hypothesis.
K-means is perfect here. The issue with K-means is designating the number of groups, but here thats irrelevant because its not the question you are asking. You are asking homogeneity and if your hypothesis is true no matter how many groups you designate for K-means, e.g. 2,3,4,5,6,7,8,9,10 (or maybe a more reduced sampling) homogeneity between target (treatment?) and general population is observed. This would therefore be a definitive analysing ASSUMING that you don't encounter a group where there is exclusively one population and the other is absent. Obviously the number of K-mean groups can't exceed the number of treated/target patients.
It would be a really good application of k-means because it avoids the question,
So why was that a priori selected?
Personally I would first present your preferred a priori group size and then say same result for all group sizes. Its simple, recognised and neat and easy for a reviewer to follow. There are ML solutions to K-means, but I think "keep it simple".

This makes allot of sense, and just knowing to google homogeneity has led me down what I think is the correct path, I am now struggling to find R packages that let me do this with explanations on how to sue them. I think clValid will work

There is a direct function called 'kmeans' in R, but its not an approach I'd use unless I was reiterating within a data pipeline (R is great for automation). For this a general stats package SPSS or Stata is easier. ML k-means are in Sci-kit learn (Python), but this is really complicated and lots of things can go wrong. Glad this helped.

source: Stack Exchange

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

sort: votes


text:
Genome annotation of a Bacillus strain
I finished the sequencing of the genome of a biocontrol activity strain Bacillus amyloliquefaciens by Illumina hiseq. I have already gone through some papers but there were no details of processing in it.
Is anywhere I can find some resources to get my research going on to genome annotation, compare it with other microbial genome data, and analyze non-ribosomal peptide synthetase, find out more genomic information about biocontrol activity?
For example, how to locate 16S DNA and some lipopeptide antibiotic genes (fungycin gene, surfactin gene) in my bacterial genome, and compare it with the reference strain and draw a phylogenetic tree. And how to use my bacterial genome and 2 reference genomes to draw a genomic circle diagram.

Welcome to the site! What is the depth of your sequencing? Did you assemble it already? Is there any reference genome of your organism? Also it would help to give a better answer if you showed what have you searched/tried.

Hi, thanks for your answer. The depth of my sequencing is 100×. The raw data is 2G, and has already assembled by sequencing company. There are two reference genome which i want to compare with, Bacillus amyloliquefaciens
FZB42，Bacillus subtilis 168.

Ok, this information is useful. Now the most important thing, what have you done to compare the annotations/annotate it yourself? You mention you had some problems finding resources, could you please provide the search terms you used, so we can avoid them and provide help about how to find the resources?

You said that you want to compare with other microbial genome data. Have you already annotated your genome? If not, you might consider a pipeline I've developed, I can gladly help if you face some problems using it.
It is called bacannot, it does a complete annotation using Prokka, Kofam, CARD, VFDB, and other resources. It might help you to have a first scratch of your data since it also collapses annotation information in a genome browser

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


text:
Short reads vs long reads with regards to break points
My professor in class said the following statements.

Short reads can map break points
Long reads can disambiguate repeat related issues

I'm very new to genomics, so I don't quite understand the first statement. What is meant by short reads can map break points. What are break points to begin with?

Did you ask your prof.? He/she is quite likely to answer your questions.

The professor's explanations are not quite clear to me and only confuse me more. Hence, I figured I ask here.

Short reads can map break points

With short reads you do an "indirect observation/inference" of a structural variant: based on your paired-end reads you figure out something is wrong because their distance or orientation is not what you expect. With long reads, you often get a direct observation of a breakpoint, with sufficient long reads on both sides (either split alignment or spanning the SV) allowing almost nucleotide accuracy.

Long reads can disambiguate repeat related issues

Short reads cannot align unambiguously in many repetitive regions, because the repeat is longer than the read. Long reads, however, span most common repeats entirely, allowing you to unambiguously assign a location for each read, as such resolving the issues.

I'd say "paired short reads can approximately locate breakpoints", because short single end reads will not do a good job, since they can't detect a breakpoints unless they land right on top of it.  Paired end reads will tell you approximately where a large structural variation is, but won't give you base resolution, unless the read happens to sit right on the breakpoint.

A "break point" typically denotes one end of a structural variation, such as an inversion or deletion. You can simply think of it as the position where the normal chromosome was broken to allow the variation. One can find these with short (Illumina) reads relatively easily since with paired end sequencing you'll have pairs spanning these breakpoints with odd relative orientations or insert sizes or such. Note that this is not a particular strength of short reads, since long reads are probably even better at this (they are not as good for finding SNPs at the moment though).


text:
What is a quick way to find the reverse complement in bash
I have a DNA sequence for which I would like to quickly find the reverse complement. Is there a quick way of doing this on the bash command line using only GNU tools?

Thanks to Manu Tamminen for this solution:
echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev

And if I want to put it into a variable instead of echo it?

I would enclose the code in $() and assign to a variable assuming you are using a recent version of bash.

that works thanks!

Reverse complement FASTA/Q:
seqtk seq -r in.fa > out.fa

https://github.com/lh3/seqtk

You may either install the emboss package (http://emboss.sourceforge.net/download/) and use the revseq program.
It contains many useful command-line programs such as the pairwise alignment tool needle.


text:
How is that ensemble Free energy is lower than the MFE in RNAfold from ViennaRNA
In the tutorial of RNAfold it states,
RNAfold reads RNA sequences from stdin, calculates their minimum free energy (MFE) structure, prints the MFE structure in dot-bracket notation and its free energy to stdout.
In the example using the -p --MEA options, the

MFE = -5.00 and the
ensemble Free Energy = -5.72.

How can the free energy of the ensemble be lower than the minimum free energy?

@MichaelG. First, thank you. I read and learn from your comments and answers. Second: I tried with CUACGGCGCGGCGC and I get: structure=.....(((...))) with -0.40 kcal/mol for the MFE and for the probab. representation of the ensemble, structu=.....,,,...)), with -1.10 kcal/mol.
I don't see pseudoknots there (should I?).

@MichaelG. 1st: Thank you! I feel honored!! 2nd: I think (and I could be wrong) that the only tool in Vienna that predicts pseudoknots is RNAPKplex (I am too far from being expert in the topic, I just googled it -couldn't find that info as fast in Vienna documentation-).  3rd: I am perplexed... I don't know how is this possible


text:
What is a signal 'chunk' in the context of Nanopore sequencing?
Does one chunk correspond to the raw signal of a single read read? Why does messing with this parameter affect basecalling results?

Hi Nick, sounds like a valid question, maybe it would be helpful to link a page that talks about this parameter. Also, there are more basecallers for nanopore, no? Which one you use?

If you have access to the nanopore community website you can read more about this here: https://community.nanoporetech.com/posts/not-chunky-enough (I'm not sure to which extent I can copy paste from there so I won't).
A TLDR version: basecalling doesn't happen simultaneously for an entire read, but it is broken in segments which are processed separately and combined afterward. Those segments are chunks. Basecalling results are affected mostly by low-complexity regions, in which a heavily biased region (in terms of nucleotide content) leads to artifacts in the per-chunk basecalling. Those issues are solved by taking a chunking-size larger than most low-complexity regions.


text:
removing mitochondrial read and unassembled "random" from multiple bam files
This is the one i used for a single bam file to filter its mitochondrial as well as unassembled read
samtools view -h SRR2920466.bam | awk '{if($3!= "chrM" && $3 != "chrUn"){print $0}}' | samtools view -Shb - > this.filter.bam

Is there any way i can do it on multiple files with parallel as i do to index and sort bam files using samtools .
So i tried this ,which is of course not working as its not able to pipe properly because it shows me error regarding additional fonts
ls *.bam | parallel "samtools view -h {} | awk '{if($3 != "chrM" && $3 != "chrUn"){print $0}}'| samtools view -Shb {.}.bam"

Any help with my above codes? As i get it there is wrong with my pipe .
After updating the samtools to version 1.9 i ran the small script
i get this error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string

I have 30 bam files i tried with 2 to see if its working or not , after the suggestion of @WouterDeCoster i tried the above one after updating samtools . the basic samtool with egrep works but for multiple file i couldn;t do it.

Your awk command is inside double quotes ("") which means that all of the awk field variables (e.g. $3) are being interpreted by the shell. Since $3 is not set in the shell, your awk is trying to compare an empty string and that's why you get the error:
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error

Note how the error message shows if( != chrM and not if($3 != chrM. Since you need to use double quotes here because you need both an outer (for parallel) and an inner (for awk) quote, and single quotes can't be nested, you will have to escape the awk variables:
ls *.bam | parallel "samtools view -h {} |
    awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb {.}.bam"

looks like i need to clear my fundamental  awk concepts let me run and will update you. now i ran your code i get output like this "�&Ct�wh���x�Y#��J�f���{/.i��7��%M�������C�̛�<kbb�������B��f����<#?���
:P�y���
u��O��r.�_2ؑ�xu�zl�,xϭ"2��Uo�ȩ�@ayi���
/0���/���۶�
�>�Xh� ��Ί�5h�f~���euSa1����;5��{���R�3L����p%��[<�1)M�y!�w��#���_�`�V5rj��v�}�
S-���y)�Nμ4��Jӕ�D���n�1��ޥ�g�����y���䒋�N��)���<:���mʲ�溎[�ʘ���K��[punit@localhost chrm]$ 62;c62;c62;c62;c^C
"

@krushnachChandra well yes, of course. You are printing a bam file which is compressed binary data. What were you expecting to see? If you want sam output, remove the b from the options you give to samtools view. The b makes it print bam output.

so i idea is after filtering the mitochondrial read and unassembled one i want  bam output written to file with bam format which will be using for downstream analysis.I don't want to print it on the terminal . tried it with -o to write the output seems like i couldn;t make it work..

@krushnachChandra -o should work. Also, you could redirect in the parallel command: ls *.bam | parallel "samtools view -h {} | awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb  > {}.bam". But as I said under your question, the new error is completely unrelated. You probably have a malformed bam file. But please ask a new question about that. Even better, run for bam in *bam; do samtools view -h "$bam" | head -n1 || echo "$bam file is corrupted"; done That should help you find the bad file.

will run it and see

its working ,the output is written to bam file

i have accepted the answer sorry i missed it to mark since i saw this thread active again .


text:
How can `limma` be used to determine genes predictive of a given phenotype?
Using limma within Bioconductor here, which is a linear model package for microarray data and lmFit on gene expression data, I can perform a series of regressions gene ~ covariates by specifying design = model.matrix(~ covariates) and calling lmFit(eset, design).
How can I instead run the regressions phenotype ~ gene + covariates?

I think you can ask this question better at bioc support site. There the developers of limma can answer you.


text:
Genotyping technologies do not maintain phase?
Why exactly does genotyping not maintain phase? I guess I'm more confused on what's the difference between sequencing and genotyping, and as a result how phase is maintained in one, but not the other. This is an image trying to do phasing using genotyping methods.
I understand how phasing works, I just don't quite get what is the difference between genotyping and sequencing, and why one of them doesn't work.

Phasing means that you are able to distinguish the both copies of your target sequence. (Usually you have one from your mother and one from your father).
If you can do this, depends mainly on the sequencing technology you are using. For example in case of sanger sequencing, you cannot do this, because you (pre)amplify all your copies at once and detect the sequence at the same time at the same location.
In case of ngs you have read information at the end. One read represents the sequence of one original dna molecule. If you detect multiple variants in one read, you can be sure that they are from the same chromosome and thus they are in phase. Now you can go on and try to find other reads that have unambiguous overlaps with these variants on your read. If you find some, you can also say that these variants have their origin on the same chromosome.

This doesn't quite answer my question. I've editted my question ot make it more clear.

Many (small) variant calling algorithms use pileups to genotype variants position-by-position. This approach alone is insufficient to phase adjacent variants. However, if the variants of interest are close enough to be covered by a single read, it's certainly possible to phase them even with <100bp Illumina reads. In fact, I was implementing just such an approach recently when I asked this question about accessing bases by aligned position rather than by read position.

Genotyping (there are of course many different methods) is likely to give you information in the format of
Locus 1: 100% A
Locus 2: 50% C, 50% T
Locus 3: 50% C, 50% A
That's not phased.
Short or long read technology will give you data in the form of reads, each of which originates from one molecule, and if a read covers two loci, you'll know how they are phased.  You will know if one DNA molecule has ACC, or ACA.
Sanger sequencing will give you mixed peaks at the locus of a mixed SNP, so there will be no phasing data.


text:
Testing for viral/bacterial infection FASTQ files
What's the best way to test if a pair of Illumina FASTQ files from blood DNA extraction contain any reads are from viral or bacterial infection from human samples? Thanks in advance.
EDIT:
Looking at the available options, I started testing Kraken, the metagenomics classification (Wood and Salzberg Gen Biol 2014) for samples where I have seen an indication of viral reads. I will update with results / further questions.

There are some people voting to close this question as too broad. Perhaps you could indicate what are your requirements and what options do you consider

If you're talking RNA, I have a blog post from a while back about Diamond/MEGAN:
http://achri.blogspot.com/2017/04/data-for-free-viral-metagenomics-in.html
In principle this will work for DNA if you have a non-trivial number of pathogen reads (since they tend to be gene dense in their genomes). Otherwise I use Kraken on a database for all the NCBI nt database, but this is slower and doesn't let you sanity check the false positive taxonomy assignments as well as I describe in the post. Of course, ideally, you run multiple methods if you have the resources to do so.

I think you could try pathoscope for this.

To know if there are bacterial genomes in the sequenced sample I have used an approach based on kmers analysis. An useful tool to do that is kmerfinder which uses a kmers database built from bacterial genomes. Quoting from the kmerfinder web site:

The method is used to find the best match (species identification) to
  the reads in one or more fastq files or one fasta file in a (kmer)
  database produced using the KMA program.

I haven't use it for viruses genome analysis, but conceptually it should be work as well.
Hope this would be helpful.

You may find this blog post helpful: Identifying viral sequences in TCGA data using Kraken and Centrifuge.

Our goal was to evaluate the performance of Centrifuge and Kraken for
  detecting viral transcripts across a cancer cohort. Centrifuge and
  Kraken produce very similar results and consistently detect HPV
  strains in cervical cancer.

There is also a review Overview of Virus Metagenomic Classification Methods and Their Biological Applications, but they focus specifically on viruses:

We evaluated 49 published computational classification workflows for
  virus metagenomics in a literature review. To this end, we described
  the methods of existing workflows by breaking them up into five
  general steps and assessed their ease-of-use and validation
  experiments. Performance scores of previous benchmarks were summarized
  and correlations between methods and performance were investigated. We
  indicate the potential suitability of the different workflows for (1)
  time-constrained diagnostics, (2) surveillance and outbreak source
  tracing, (3) detection of remote homologies (discovery), and (4)
  biodiversity studies.

What type of NGS data is it? If smallRNA or RNAseq data, you can simply assemble the reads into longer contigs and do the blast search using your reference viral or bacterial sequences.


text:
PLINK - Transposed BED file?
I'm currently working with the PLINK file system to store large amounts of genotype data.
The plain format consists of three files, two files for storing phenotype and marker information as well as one file (.PED) containing the genotype matrix.
There exists a binary version of the .PED file (.BED), which stores the data in a number of consecutive byte blocks, with each block encoding the genotype of all samples for a specific marker.
While there exists a transposed version of the .PED file (.TPED) which is also supported by PLINK, I couldn't find out whether there also exists a transposed version of the binary .BED file.
Does it exist at all? Are there other file formats to compactly store a genotype matrix for fast sample / marker extraction?

Yes, this exists, and can be efficiently generated by plink 2.0's --export ind-major-bed command.  (The third byte is 0 instead of 1, and the specification is otherwise identical to that of regular plink .bed files with samples and variants swapped.)

Thanks you. So the sample/variant-swapped binary file has the same .BED suffix? Does PLINK2 detect the file based on the third byte and is there a noticeable speedup when extracting a samples genotypes?

It has the same suffix, yes.
PLINK 1.9 and 2.0 automatically create temporary variant-major .bed files when they see a sample-major .bed file, since almost all the code assumes variant-major, so this never speeds up PLINK.  But it can be useful input for other programs.


text:
Knowing which positions are variants
Let's say I'm sequencing a donor, and I have 20X coverage. And let's say all the reads covering a particular base has the same allele at the position. Then, is that position considered a "variant" with only 1 allele in that position or is that not even considered a variant?
The reason I'm asking is, I was given some read data, and it's in the form location:variant - location:variant.
But there was often the case where for all the variants at a particular location, there was only one allele. But this data is labelled to only have variants. Hence, how is a particular location deemed a "variant."
For example, lets say I'm told that the following is "variant reads given in the form read# location:variant - location:variant"
read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

Then, why are these two locations considered variants? It seems like it's just a heterozygous locations. Why were these locations considered to carry variants? Similarly,
read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

In those reads, it just seems like it's homozygous. Why is it being called a variant? How was it classified as such?
Put another way, if I'm given 100000 reads, and I said I only want to look at informative reads, and so I'm given a subset with only 100 reads, and those 100 reads look something like the above, then how were these reads selected?

Please [edit] your question and include an example. I am afraid it is very hard to understand what you mean.  I suspect you are using the terms variant and allele in a non-standard way, so an example illustrating what you mean exactly would be really helpful.

@terdon I've editted my question.

Whether you call something a variant or not depends on what you are using as a reference. You need the genomic position, not the position with respect to the read and you need a reference genome to treat as a "truth set". You cannot know whether something is a variant or not just by looking at the reads without a reference genome. What kind of data are you analyzing? Have you aligned your reads to a genome?

I see. These position, I know, are positions on the genome. Hence, I'm guessing these reads have been mapped to the reference, and then being called variant. Is that correct?

So, in the 2nd read cloud, if the reference was heterozygous, but the donor reads are telling me it's homozygous then it's called a variant. Am I understanding that correctly?

Essentially. It's a variant if what you have in your sample isn't what you find in the reference at that position. So, if the reference is A/A and you have A/T, T/T or T/A then it's a variant.

Set 1
read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

I understand this as showing the location of 2 variants, supported by 2 reads each. Now, if the reference has an A at position 6523432, your sample could be heterozygous for an A->C variant at position 6523432. If the reference has a C at that position, then your sample could be heterozygous for C->A. If the reference has something else again, for example G, then your sample could be heterozygous for G->A/C. Similarly for location 6710432.
Set 2
read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

In this case, if the reference has an A at position 6523432, the reads do not support a variant. However, if the reference has anything else, then the reads would support a homozygous variant e.g. G->A. And the same for position 6710432.

Now, in all of these cases, you cannot know whether a variant is real without taking into account the number of reads that support it, the quality of the alignment, and the quality of the reads. That's why we do this using dedicated variant callers, software that is designed to take all these factors into account and decide whether to call a variant at that position.
But in all cases, the variant must be defined with respect to a reference genome. You cannot call a variant just by looking at reads without any knowledge of a reference. A "variant" just means "not what is found in the reference". It is largely a matter of random chance whether something is a "variant" or "normal". If the individual(s) whose genome was taken as a reference has the same sequence, then you call it normal. If not, then you call it a variant. There are many cases where the standard hg19 or hg38 genomes have a "normal" residue which is actually far less common than the known "variant". There's nothing special about the reference genome, after all, it's just the one we chose to call a reference.

Perfect, this exactly the clarification I needed, thank you!


text:
What could cause differing counts of R1 and R2 in Paired End Sequencing (RNASEQ)
I recently finished mapping an RNAseq run using STAR2.3.0 and noticed the read1 and read2 counts are different, according to samtools flagstat.
The map% is ~80% but the R1 R2 counts are:
R1=139283692
R2=137472495
This usually does not happen for PE sequencing, which makes me a bit nervous to proceed with analysis. What are some potential reasons something like this could occur? I do not have experience with this happening and can't seem to find reasons by searching the web.

Presumably, the actual read counts in your two fastq source files are the same, right? I have no experience with RNAseq but the first thing that comes to mind is a problem in the sequencing step where too many reads were of low quality and were discarded by the aligner. Does that make sense?

The actual read counts of the two fastq source files are the same. It makes sense and was a thought of mine. I was wondering if there is a concrete way to check this out?

Use FastQC to compare the overall quality of both fastq files.

@d_kennetz Does the actual read counts matches R1? Have you checked which of the reads are missing in the bam file?


text:
Proteins with one SS bond?
I would like to find proteins with exactly one SS bond. Is there a database where I can search this? I've tried advanced search on https://www.rcsb.org/, but no such option, at least I could not find it.

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.

Thanks! I need just a few examples, to the 1st method is perfect!


text:
How do I find split reads?
How can I detect a split read in a BAM file?
Is there any sign in the CIGAR string that describes split read?

Welcome to the forum.  Have to done any research on this topic?  Have you looked into secondary and/or supplementary alignments?  What aligner was used to generated the BAM file?

yes, for example this site " https://galaxyproject.org/tutorials/ngs/" describes SAM/BAM file structer but i want know what is "split read" and how reported in this file?

@fatemehkhoramdel That answers (to an extent) the question on "have you done any research?". Can you also answer the questions on secondary/supplementary alignments and the tool used to generate the BAM please?

If i understand this question correctly, you can look at  "The Sequence Alignment/Map format and SAMtools" at google scholar. otherwais i don't know about it

It's easier to help you if you provide a small amount of background info in your question to see what you know and what work you've done so far. Another good first thing to do is to search the Bioinformatics StackExchange to see if a similar question has been asked before. In your case, it looks like this question may be useful for you

Hello @fatemehkhoramdel, I find this a nice question. The context you write in comments is very useful, it would be even better if you would [edit] your question and it the info there directly :-).

The answer to this question depends on the aligner used to generate SAM/BAM file. For BWA method that is part of the pipeline mentioned in comments, split reads will be flagged as supplementary alignment (2048). You can modify this behavior by using option -M as described in the BWA-MEM manual page.
So, in case of split alignment, you will end up with multiple aligned segments that have identical query name and they are flagged as 2048. These segments may only have part of the original read sequence, but the complete read sequence is present in at least one of the segments.
Also, there is no unique sign in CIGAR of an individual segment to identify split reads. Each aligned segment of the original read may have long deletion in its CIGAR or aligned pairs. However, without looking at all segments with the read ID it is not clear whether the segment has a novel sequence in it that doesn't align with reference or that part matches different chromosome/coordinates.


text:
Building cancer gene database
I am trying to build a database of different cancer genes i.e lung cancer, throat cancer etc. Searched through NCBI website by searching for nucleotide sequences, but I am unable to find appropriate results. Any help in this regard is much appreciated

COSMIC has a set of mutations and the Cancer Gene Census already curates cancer-related genes. What is the purpose behind your venture?

@RamRS To perform clustering

Why not just cluster based on existing annotations? You could get information from the CGC, use the HGNC symbols or other identifiers to download additional data using biomaRt and then perform all sorts of analysis on that. Why build a new database?

@RamRS i have taken up this as part of my PhD , wherein I will be clustering cancer gene data treating them as pure categorical data and ultimately try to identify similar patterns across clusters if they exist

What will you be clustering (the term data is too generic and vague) and what sort of patterns will you be looking for? If this is part of your PhD, you definitely would have looked into literature on such protocols and the sorts of data available at the gene level on various cancer types. Can you give us any more detail on it?

@RamRS gene expression to be specific

And will you be clustering gene expression data from a public resource such as TCGA? I ask because if it's based on public datasets, this would not be a novel endeavor at all. In fact, there exist explorer tools that let you run correlations between various omics dimensions on the fly.


text:
Find similar proteins with similar tertiary structure?
I'm a mathematician studying knotted curves and am trying to find some applications of my work in proteins. Using https://www.rcsb.org/, how do I find proteins that have similar  (but perhaps not the same) tertiary structure? Preferably they would have similar primary structure and function...
Perhaps what I'm looking for is "similar proteins", if there is such a thing?
Also, is there a way in the database to select more "common proteins" (important or well studied, not just some rare exotic one (if there is such a thing)?

As a general question on the board, more detail is often good, because of the answers can be inaccurate due to misunderstanding of the biological context.
There are several ways to do this and at least 3 that I'm aware of.

Blast - sequence similarity searching the PDB database
Contact map identity, I use this approach, but its complicated - and its precision or utility is context specific, i.e. it could be loads of work and miss the vital detail you need.
Direct mapping one protein structure against another, this in particular is true for specific domains, or sites within the protein and the result is angstroms. Its not difficult and I assume this can be scaled to the whole protein.

Point 1 "Blackboxes"
By far the easiest is option one and I suspect I-Tasser (here) use this principle before instigating "threading". This package is considered "good" (others appear to consider it the best of its type), I believe it has done well in CASP tests (Critical Assessment of Techniques for Protein Structure Prediction). So the simple answer is dump your sequence in there an watch the output. Pyre is something similar. These are interesting tools but certainly not definitive, really there a starting point.
Point 1 "Formal Blasting PDB"
This is the easiest and most mature approach for your question, without knowing the details. I've run a test to ask the following question.
Which protein structure available on PBD is most similar to Japanese encephalitis virus (JEV) based on protein homology?

Step A
I've no specific interest in the virus BTW. The most robust Blast engine IMO is NCBI here . I dump a JEV sequence into the box, and click three important options (otherwise it doesn't work).

Database = protein databank (PDB), i.e. structural data
Organism "Japanese encephalitis virus (taxid:11072)": the taxid will come up automatically (if not just proceed as normal)
Next to the Organism box I click "Exclude", otherwise this whole approach doesn't work
Obviously I select "blastp", otherwise it definitely will not work

The output is/was here (but is now inactive). Summarising what it presented, E is so high its "0" (>2e-139) protein ID "77.59%" and just states "envelope protein" Chain A. I hit the link and the protein that comes up has a name which ends in "6A0P_A". This is the PDB reference and Percentage ID and PDB reference are the only two bits of info you need. The PDB reference needs the "_A" removing, thus its 6A0P. The _A just refers to chain A of a dimer.
Step B
You then go back to the Web site you mentioned and input the PDB ID you stated into the search bar. This results in this quite interesting page here. The structure is Usutu virus, which surprised me.
Point 2 . This is much harder and not only involves specific Python code, but a post unsupervised learning analysis and results in an average similarity of structure over the whole protein. There are a number of published algorithms for generating a contact map, or its easy to write your own. Its more advanced, its complicated and unless you've seen loads of them they not immediately understandable. Contact maps angle at a very specialist approach of predicting and understanding protein structure. It is not for the faint hearted and some biological details will be omitted, which may or may not be important for your purpose.
Point 3 . I don't have research grade expertise here, so I forget the exact button to click, but its quite easy via PyMol and other graphical viewers to perform spatial similarity for specific epitopes and this is something you need to learn within your sphere. Basically you will load two proteins into the image viewer and then they will be a "similarity aligning" function and you zoom in on your output to generate an answer in angstroms. You will need to learn this approach for the later parts of your project.
Summary, The orthodox way would be to perform point 1, select several targets and then examine them in detail via point 3. I definitely recommend direct Blast if you go down this route, because it is the most tractable and avoids potential unknown subsampling.  I personally would not recommend Point 2 in a general context because you have to consider an unsupervised learning approach towards processing a heap of contact maps. Ultimately, which ever approach you take you can never avoid Point 3 (which someone else does, in my case).
Note This may or may not be obvious to your but just because two sequence have the closest homology doesn't mean they will have the closest structural similarity.

An explicit model of a "knotted curve" would be very interesting by the way because it would enable an analysis of hinge regions, which are quite important in viruses. I can immediately think of a entire application thereafter via an algebraic solution based around modelling complex non-linear lines.


text:
scaling algorithm of deeptools
I am using deeptools to generate a metagene plot. I am wondering how the deeptools scale genes if the size of genes is shorter than bin size. For example, if my bin size is 1000 bp and the gene is about 450bp, how does the deeptools computeMatrix scale it to make a metagene plot of 1000 bp?


text:
Delly finds identical deletions with different consensus sequences
I am using delly to call structural variants on a human WGS sample. I noticed that it finds certain variants that are essentially duplicates. For example:
chr4    177594033       DEL00005339     G       <DEL>   .       PASS    PRECISE;SVTYPE=DEL;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chr4;END=177594095;PE=0;MAPQ=0;CT=3to5;CIPOS=-3,3;CIEND=-3,3;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=5;SRQ=0.993243;CONSENSUS=GACAAAAACATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAAC;CE=1.75466     GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-51.3797,0,-42.791:10000:PASS:1115:1821:1820:1:0:0:23:27
chr4    177594033       DEL00005340     G       <DEL>   .       PASS    PRECISE;SVTYPE=DEL;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chr4;END=177594095;PE=0;MAPQ=0;CT=3to5;CIPOS=-3,3;CIEND=-3,3;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=12;SRQ=0.994186;CONSENSUS=ATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAACAATGTAATAAGTATGTTTGAATAGTATTATGAT;CE=1.72637    GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-51.3797,0,-42.791:10000:PASS:1115:1821:1820:1:0:0:23:27

The two variants above have the same start position and the same end position (END=177594095). However, they differ in the SR, SRQ and CONSENSUS fields. These are defined in the VCF header as:
##INFO=<ID=SR,Number=1,Type=Integer,Description="Split-read support">
##INFO=<ID=SRQ,Number=1,Type=Float,Description="Split-read consensus alignment quality">
##INFO=<ID=CONSENSUS,Number=1,Type=String,Description="Split-read consensus sequence">

Of course, the consensus sequences are quite similar, but they aren't identical (I have simply added some spaces to the beginning of the 2nd sequence here to align them):
GACAAAAACATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAAC
         ATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAACAATGTAATAAGTATGTTTGAATAGTATTATGAT

So what am I looking at here, exactly? If this is a deletion on chromosome 4 from position 177594033 to position 177594095, then how can the length of the consensus sequence be different? This is the "split read consensus", so presumably the sequence that appears to have been deleted based on how the reads covering that region are split. This means that these reads align to the genome with a gap and this gap indicates the deletion (image source):

Therefore, if the start and end position are the same, the consensus should also be identical. And, in any case, if the start and end positions are the same, then this should be the same structural variant, so why is it being reported twice? How can I interpret this result?

It turns out that this was just a bug which has been fixed by this commit. The current master version of delly no longer exhibits this behavior.
If you are facing this issue, either build the current master version or wait a bit for the next release which should also include this fix.


text:
Getting this information for GISTIC2
Long time I am struggling how to provide GISTIC2 required input like Num_Probes and Segment_Mean
I have called somatic copy number by scatngs and I have copy number output file (attached);
https://www.dropbox.com/s/r6nef2v1j2di6q6/my.copynumber.caveman.csv?dl=0
This is few lines of scatngs copy number output
1,1,13116,665401,2,1,3,1
2,1,668374,96324627,2,1,2,1
3,1,96326095,96419709,2,1,1,0

In scatngs reference paper says in this file we have
Segment number

Chromosome

Start position (origin-1)

End position (origin-1)

Major copy number—normal

Minor copy number—normal

Major copy number—tumor

Minor copy number—tumor

scatngs also returns another output like below
https://www.dropbox.com/s/7i28vooixdp06c2/WTSI-OESO_121_1pre.copynumber.txt?dl=0
Chromosome  Position    Log R   segmented LogR  BAF segmented BAF   Copy number Minor allele    Raw copy number
rs62635286  1   13116   -0.0290741822638953 0.245246748556562   0.25    NA  3   1   2.35806652809599
rs75454623  1   14930   0.402869048014209   0.245246748556562   0.3095  NA  3   1   4.14031506947036

I am wondering how I can get Num_Probes and Segment_Mean From this files?
People say
Segment_Mean = log2(tumour copy number / 2)

So, does that mean I should add up Major copy number—tumor and
Minor copy number—tumor to get tumour copy number? for getting Num_Probes what should I do?
EDITED

people says
scatngs uses read counts at SNPs to generate calls and so the
  locations of the SNPs is what I need to
  generate the number of markers in each segment. (eg each SNP = 1
  marker).
Do they mean, for instance if the position of SNP  rs62635286 is
  13116, so the num_prob is 13116?

Any help please?
Thanks a lot in advance

Cross-posted on GitHub for some reason: https://github.com/cancerit/ascatNgs/issues/84

I mentioned a cross-post. That does not stop anyone from reading the question or helping you.

Mentioning a cross post is hardly negative feedback or reason to ignore a post that belongs on bioinfo SE.

Did you try reading the GISTIC2 documentation about how to get the data in that format? Or perhaps they point to a known file extension (or a weird one from only a single program) I tried searching on a couple of website about these required Num_probes but couldn't find anything. Are you using this software or other with the same name?

I read GISTIC2 documentation, there is nothing about how to get required format. Exactly I am trying the same software you are referring to. People using GISTIC in broad institute produce their CN by this method that returns format for GISTIC https://software.broadinstitute.org/gatk/documentation/article?id=9143#1

Depends on the software to generate the segmentation file;
For example I have used scatngs for segmentation that for sample returns
Segment number
Chromosome
Start position (origin-1)
End position (origin-1)
Major copy number—normal
Minor copy number—normal
Major copy number—tumor
Minor copy number—tumor

By  log2(Major copy number—tumor +  Minor copy number—tumor /ploidy of each sample) you would get segment mean
For number of probes
The number of probes depends on where the data came from.
1) From SNP arrays:  Number of array probes within each segment.
2) From exon capture data:   Probes are the ends of exons and so the number of probes are equal to the number of exons  X 2.
3) From whole genome sequencing:  Probes are the boundaries of bins used to generate the copy number data and so the number of probes are equal to the number of bins +

For example in my case for scatngs
Each SNP is a marker, and so you have to count how many of these SNPs are on each segment

I have asked this question in
Extracting the number SNP in each range


text:
How to compute LD for pairs of variants with Plink
I have a set of SNPs (single nucleotide polymorphisms) { S1, S2, ..., SN } from approximately 200 humans. I wish to determine the linkage disequilibrium (LD) for a defined subset of these SNPs. By this I mean e.g. I would like to compute e.g. LD between S1 and S2, S3, S4, but not the LD between S2 and S3.
Could this calculation be performed in PLINK?

Hi Learner, thanks for you post and hope you can hang around here.The obvious question is why do you wish to sub-sample the genome, i.e. biological rationale?

In this case, I have a priori generated SNP pairs as the output of an eQTL analysis. Now I just want to know about the LD between pairs of SNPs, which I could compute manually, but I have thousands of pairs, and so it is ideal if there is efficient software like PLINK that I could use...

I've added a comment below

It could be that your question is not clear, but reads that you are wanting some SNPs subject to LD, but others excluded. For an LD to work you'll need a large number of SNPs which would be input simultaneously. You can subsample a given (human) population, your data set might therefore be further broken down from 200 humans.

Yes, you can do this with plink.
You can use the flag --ld-snp <variant_id> to calculate LD between your SNP of interest and all other SNPs.
If you have multiple SNPs of interest, and want each of those compared with every SNP in your dataset, use --ld-snps <snp1-snp100,snp102,snp104>
A final option to calculate LD between a specific pair of SNPs is --ld <variant_id1> <variant_id2>
This the relevant plink manual page

Anyway, it would be fairly easy via the filtering options here
The --exclude-snp  is what you are wanting.

I have worked in LD and subsampling a genome is open to forcing a result. The alternative approach is the subsample is the population, or examine disequilibrium with respect to the physical chromosome. A local shift in disequilibrium  or heterozygosity  (sliding window along a chromosome) is an indication of selective pressure. This approach is frequently used in eukaryotic parasites to identify drug resistance. I don't have experience of human genetics but there are certainly enough here who do. However, I think you need to define the biological rationale for the subsample. Just having loads of data is easily resolved because any algorithm can be written to handle automation.

The issue is that LD usually involves a 5% cut-off, even using a maximum likelihood approach I'd be cautious (I don't which PLINK uses), and is therefore sensitive variations in the data set


text:
Question about enrichment score (ES) and NES from GSEA
I have two results from two separate GSEA runs (2 different inputs / samples with same gene sets), and i want to know which sample has the most enrichment about specific pathway (gene set).
I want to compare the ES from the two different runs, looking at the same gene set. Can I use ES (or NES) from two different GSEA analysis about same gene set?
or should I combine the input dataset into one, and run the analysis again?

Usually GSEA is not applied to samples but to comparisons. Could you further explain what is your input and the method used for GSEA? Perhaps it is one I haven't heard of (there are many enrichment methods not only the one from the broad institute)


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

What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?
I saw some papers are doing this, is it to validate the docking protocol?
Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

Usually the Km I've seen related to K_cat is the michaelis-mendel constant, not a docking score (but maybe this is different terminology). Could you explain some context of this question? A text of the paper or a reference where do you find this. Also if you could show what you have searched/read it would allow us to skip these resources when trying to find the answer.

https://www.nature.com/articles/ncomms12965

refer to the graph in page 3

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

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

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

Totally. Docking finds were ligands bind. Although the energetics don't match too well the experimental kinetics as of 2019.
The catalytic efficiency (k_cat/K_M) is the ratio of the turnover (k_cat) over Michaelis constant (K_M). The Michealis constant is in turn the ratio of on-rate over off-rate plus turnover. The ΔΔGibbs of the ligand ought to correlate with the ratio of the on-rate/off-rate, while the turnover is unrelated.
The turnover should correlate with the difference in energy between the substrate bound form and the transition state. This is the realm of theozyme design and requires correctly protonated ligands and protein, but when it is you get a nice transition state hill (e.g. personal example).
There are a few "ought to" in the above. Docking and force-field energy scores are not precise. The paper you refer to is doing very advance docking (akin to dynamic undocking) and, not accidentally, the ligands are very flexible and rather soluble.

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

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


text:
Why are my kallisto and salmon results differing so much just for lncRNA transcripts?
I am running some analysis on an RNA-seq dataset. I have a list of transcripts that are potential lncRNA for which I ran both Kallisto and Salmon aligners. The input data for index building and quantification includes mRNA as well as these potential lncRNA.
Commands for one replicate:
Salmon
salmon index --keepDuplicates -t transcripts.fasta -i index_directory/ -k 31
salmon quant --seqBias --gcBias -i index_directory/ --libType ISR --auxDir index_directory/ -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz -o output_name

Kallisto
kallisto index -i kallisto_index.idx transcripts.fasta
kallisto quant -i kallisto_index.idx -o output_name --rf-stranded sample_R1.fastq.gz sample_R2.fastq.gz

Then I used the tpm values, averaged replicates (3) and plotted the number of transcripts per condition that had a tpm > 2. Same result for tpm > 5.
The following graph shows the for mRNA:

And for lncRNA:

Again, the analysis were done with both lncRNA and mRNA transcripts together in the same fasta file.
The overlap is not shown here but is always smaller than the count for Kallisto, for the lncRNA.
In literature both tools seem to have similar results. I cannot find a reason for this. If there would be an error somewhere I would expect to see also differences in the conding RNAs which does not happen.

While the question seems great I doubt we can provide an answer without knowing the data. Could you reproduce this in a public dataset? (I think you won't be able to share the (processed) data isn't it?) But perhaps the errors lies somewhere in the code used to create the plots (rare, but it could happen)...

It is true, I am not able to share details on the data. However I can try to reproduce the results with a public dataset. Not just to edit the question, but it can actually give some clues on what is happening in case of a code error or due to the type of data. Thank you for the remark. I wanted to see if anyone had come across this problem since it is not the main focus of my project so I am not sure when I can adress this

I think a scatter plot of the abundance estimates from each method might be a more informative comparison. That is, for each lncRNA, plot kallisto's estimate vs Salmon's.  Are the kallisto estimates systematically lower or are they non-existent?

Shh, don’t tell Lior Pachter.

There has been a detailed comparison of kallisto and salmon for lncRNA transcripts by Zheng et al https://www.biorxiv.org/content/biorxiv/early/2018/01/09/241869.full.pdf
It shows them to be near identical. This suggests there is some problem with the way you made the comparison/plot.

I have had similar results to the preprint when I did some comparisons between the two methods.


text:
How do you download stuff from NCBI fast?
I have 2 .csv files that contains a list of accession codes. For example for this experiment a .csv file will contain:
SAMN10870244
SAMN10870243
SAMN10870242
SAMN10870241
SAMN10870240
SAMN10870239
SAMN10870238
SAMN10870237
SAMN10870236
SAMN10870235
SAMN10870234
SAMN10870233

Then I run the following
for i in $(cat *.csv); do echo ${i}; /usr/bin/sratoolkit/bin/fastq-dump --gzip ${i} & done

Some of my downloads takes a lot of time days. Is there a faster way to get stuff from NCBI?
Are there mirrors of NCBI to accelerate the process?
Why is it so slow?

Answers to this question can be found on biostars. To give you a short summary:

Faster way might be to use the parallel-fastq-dump, as suggested in this answer. I never tested that tool though, my own experience is that prefetch is more stable than fastq-dump command.
In that same post it seems that the Japanese mirror seems faster, don't know if that is still valid.
The reason this is never optimized might be that downloading from SRA is not a task you'll have to do repeatedly, as described here in this post.

Thanks, comprehensive answer. Is it this site https://www.ddbj.nig.ac.jp/  ?

https://www.ddbj.nig.ac.jp/dra/index-e.html, yes it is the same.

You can try wget to download SRA files from NCBI server. This seems to be 2x faster in terms of CPU time than prefetch and has the advantage of not clogging your /home/ncbi/public/sra/ directory in case you have space limitation. The tutorial for wget downloading method is available at https://github.com/johanzi/ncbi_tutorial.
Here my timing of the prefetch method:
time prefetch SRR3300113
real    3m5.609s
user    0m6.332s
sys     0m2.076s

Versus the wget method:
time wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR330/SRR3300113/SRR3300113.sra
real    3m4.696s
user    0m0.456s
sys     0m3.332s

That makes a user+sys time of 8.402s for prefetch and 3.788s for wget. Not so much difference for real time though (about a second).

Hello Johan, welcome to BIoinformatics.SE. Is there a reason to think that wget would be faster? Intuitively I would guess that it's the server that is clogged, then using a different program to query the server does not sound like something that would help. (by the way, my inability to see how it could help does not mean that it does not :-)

Hi Kamil, the wget process is about 2x faster in terms of CPU time than prefetch method based on my test:

user + sys sort of tells you how long CPUs were active during the real time of execution (i.e. for user real is what matters, you can check this). So, what I would get out of your nice test is, that the wget is a more optimized tool, but both of them are fetching the data from exactly the same place and the limiting factor is the speed of their ftp.

Yes. I would say the main advantage of wget is that you can choose where to put your downloaded SRA files while it goes in default into your home directory /home/user/ncbi/public/sra with prefetch. It can be a problem when you work on a server and you have limited disk usage for your home directory.

You can choose where to put the cache, at least as of a few versions ago, via vdb-config. I agree it's a total pain that it defaults to the home directory instead of the current one, but it's a fixable problem.


text:
Does StringTie use exon length for gene abundance or complete gene length?
I'm currently using StringTie with option -A for calculating gene abundance (I'm interested in TPM). While checking the manual, I couldn't find any reference to the record type StringTie uses to calculate this value. Does anyone know if it is total gene length or total effective (i.e. exon) gene length?


text:
What is the best argument for maximum parsimony method in phylogenetic tree construction?
Beside it (maximum parsimony) being compuationally cheap what other good arguments are there for it? Is there any model behind this principle ? Why would one expect this principle to provide right phylogeny in any situation at all ?

You cross-posted this on Bio - cross-posting is discouraged across the SE network.

Noted. I didn't now that.

The principle of parsimony states that the simplest explanation consistent with a data set should be chosen over more complex explanations... (Stewart 1993)

Sometimes independent events in different lineages cause similarities in molecular sequences (or morphology or behavior or whatever characteristics one happens to be studying). These events are not true evidence of recently shared ancestry, but it is easy to mistake them as such evidence.
The principle of parsimony asserts that we will be better off assuming the simplest case--that is, that identical characters are evidence of recent shared ancestry, since the simple explanation should be true more frequently than alternative explanations. There are cases for which the maximum parsimony model works very well, and there are cases for which maximum likelihood, distance matrix, Bayesian, or other models perform better.

Stewart CB (1993) The powers and pitfalls of parsimony. Nature 361, 603-607, doi:10.1038/361603a0.

Alot happened in phylogeny since 1993

And yet the principle of parsimony remains the same, no?

Not entirely sure what you mean, but parsimony is still parsimony and the fundamental idea of assessing a tree as an aggregate of every individual alignment position remains the same.

The OP asks why it even works at all, and my answer is that in some cases the simplest assumption is sufficient to infer the correct relationships. I don't see anything there that is compromised by the impressive advances in phylogeny inference that have occurred since.

IMO had parsimony focused on protein trees - that might have succeeded. It gets a bit stuck with DNA. A key problem was they (parsimony exponents) didn't defend their technique rigorously, but I suspect you've never come up against a parsimony reviewer. Well they defended their technique with opinion rather than analytical fact.

I think you're mistaking me for a proponent of parsimony as opposed to someone trying to answer a question. :-) Don't worry, I was just as alarmed (and amused) several years ago when Cladistics announced it would only accept parsimony trees!

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.

Parsimony was a nice idea and developed the algorithms and tree methods towards its successors, maximum likelihood and Bayesian. The fundamental assumption of mutational independence at each alignment position was founded via parsimony. Just the recursion algorithm involved in  reading a tree isn't trivial coding. Also Walter Fitch (Fitch Marigold method) was a nice guy.
Weighted parsimony has performed well in simulation studies to prevent 'long branch attraction' artefacts. The key issue that was never resolved was how to objectively weight an MP analysis. It could be resurrected however, because machine learning ML could independently provide the weights in specific contexts. That a lot of "coulds" though.


text:
Repeatmasker output table
I'm working with an unassembled genome of cephalopod and I would to find transposable elements. My idea is to use repeatscout to generate a library and repeatmasker to find the elements. After creating a library with repeatscout,  I run Repeatmasker but it generates a table like this:
 file name: OB_100DEC.fa
sequences:           100
total length:  209466439 bp  (184235452 bp excl N/X-runs)
GC level:         35.24 %
bases masked:   73100024 bp ( 34.90 %)
==================================================
               number of      length   percentage
               elements*    occupied  of sequence
--------------------------------------------------
SINEs:                0            0 bp    0.00 %
      ALUs            0            0 bp    0.00 %
      MIRs            0            0 bp    0.00 %

LINEs:                0            0 bp    0.00 %
      LINE1           0            0 bp    0.00 %
      LINE2           0            0 bp    0.00 %
      L3/CR1          0            0 bp    0.00 %

LTR elements:         0            0 bp    0.00 %
      ERVL            0            0 bp    0.00 %
      ERVL-MaLRs      0            0 bp    0.00 %
      ERV_classI      0            0 bp    0.00 %
      ERV_classII     0            0 bp    0.00 %

DNA elements:         0            0 bp    0.00 %
     hAT-Charlie      0            0 bp    0.00 %
     TcMar-Tigger     0            0 bp    0.00 %

Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %

Small RNA:            0            0 bp    0.00 %

Satellites:           0            0 bp    0.00 %
Simple repeats:  290711     15358198 bp    7.33 %
Low complexity:       0            0 bp    0.00 %
==================================================

* most repeats fragmented by insertions or deletions
  have been counted as one element

The query species was assumed to be homo
RepeatMasker Combined Database: Dfam_Consensus-20181026

run with rmblastn version 2.6.0+
The query was compared to unclassified sequences in ".../OB_100DEC_repeats_filtered1.fa"

This was my command to run repeatmasker:
./RepeatMasker -s -lib /home/RepeatScout-1.0.5/OB_100DEC_repeats_filtered1.fa /home/Workdirectory/OB_100DEC.fa

The file OB_100DEC_repeats_filtered1.fa is the output of repeatscout and looks like this:
>R=3 (RR=4.  TRF=0.000 NSEG=0.000)
TAAGGCGGCGAGCTGGCAGAATCGTTAGCACGCCGGGCGAAATGCTTAGCGGTATTTCGTCTGTCTTTACGTTCTGAGTT
CAAATTCCGCCGAGGTCGACTTTGCCTTTCATCCTTTCGGGGTCGATAAAATAAGTACCAGTTGAGCACTGGGGTCGATG
TAATCGACTTACCCCCTCCCCCAAAATTTCTGGCCTTGTGCCTATATTAGAAACGATTATT
>R=4 (RR=5.  TRF=0.122 NSEG=0.226)
ACACACACACACACACACACACACACATATATATATATATACATATATACGACGGGCTTCTTTCAGTTTCCGTCTACCAA
ATCCACTCACAAGGCTTTGGTCGGCCCGAGGCTATAGTAGAAGACACTTGCCCAAGGTGCCACGCAGTGGGACTGAACCC
GGAACCATGTGGTTGGTAAGCAAGCTACTTACCACACAGCCACTCCTGCGCCTATATATAT
>R=6 (RR=7.  TRF=0.134 NSEG=0.247)
TTGTTTCAGTCATTTGACTGCGGCCATGCTGGAGCACCGCCTTTAGTCGAGCAAATCGACCCCAGGACTTATTCTTTGTA
AGCCTAGTACTTATTCTATCGGTCTCTTTTGCCGAACCGCTAAGTTACGGGGACGTAAACACACCAGCATCGGTTGTCAA
GCGATGTTGGGGGGACAAACACAGACACACAAACACACACACACACATACATATATATATATATATATATA
>R=8 (RR=9.  TRF=0.391 NSEG=0.000)`

and it is composed by all the TE (not classified) find in all the scaffolds of the genome.
Can anyone explain why there are no TE?
The fasta file that I used is only a part of the whole genome, it is about 200 mb and it is composed by the 100 greatest scaffolds of my genome.
Maybe I made a mistake with the use of libraries. Anyway I installed RepeatMasker with rmblastn version 2.6.0+ and RepeatMasker Combined Database: Dfam_Consensus-20181026.
thank you

Please [edit] your question and clarify. What species are you working with? Also, you mention contigs, so should we assume the genome isn't assembled? How long are your contigs? Are they long enough to enable finding LINES? What exactly do you have in OB_100DEC_repeats_filtered1.fa?

I modified my post whit all the information.

As far as I can tell, you use the -lib option for species where there is no existing RepeatMaster database. Therefore, all it does is attempt to locate versions of the repetitive elements from the -lib file in the target genome file. However, since your lib file doesn't have any information on what type of element each of your hits is, it looks like RepeatMasker is treating all of them as "Unclassified".
Note the following from the repeatmasker.help file you should have downloaded along with RepeatMasker:

The recommended format for IDs in a custom library is:
>repeatname#class/subclass
or simply
>repeatname#class
In this format, the data will be processed (overlapping repeats are
  merged etc), alternative output (.ace or .gff) can be created and an
  overview .tbl file will be created. Classes that will be displayed in
  the .tbl file are 'SINE', 'LINE', 'LTR', 'DNA', 'Satellite', anything
  with 'RNA' in it, 'Simple_repeat', and 'Other' or 'Unknown' (the
  latter defaults when class is missing).

This makes sense. RepeatMasker cannot know if an arbitrary sequence you give it is a SINE or a LINE or anything else. It will therefore rely on the categories present in the fasta headers. Since your file has no categories, it identifies everything in your lib file as "Unknown", which is presumably what the "Unclassified" section in your output refers to:
Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %

Ok, I understand the problem. In your opinion, is it possible to run with repeatmasker OB_100DEC_repeats_filtered1.fa instead of the genome.fa as fasta file and then observe the results? Could I observe any results?

@jonny no, that would just find the same thing. I would instead try running Repeatmasker directly, setting the species to cephalopods:  ./RepeatMasker -s -species cephalopoda /home/Workdirectory/OB_100DEC.fa


text:
Authoritative source on human cytogenetic regions?
I am looking for a database that would keep track of  human cytogenetic regions and genomic coordinates per genome assembly. I had expected the Genome Reference Consortium to have it, or Ensembl, or NCBI, but I could not find it.
Is there a place with a translation table between human cytogenetic regions and genomic coordinates?
EDIT
After benn's answer I found that you can also get this information from Ensembl's REST API: https://rest.ensembl.org/info/assembly/homo_sapiens?content-type=application/json&bands=1.
I have transformed the data provided by Ensembl and generated a CSV file with the cytogenetic regions:
https://github.com/ramiromagno/gwasrapidd/blob/master/data-raw/cytogenetic_bands.csv

Do you look for a database a paper with a table or what? To which genome should these genomic coordinates be from?

@llrs: I was looking for a database that would keep track of these definitions per genome assembly. I had expected the Genome Reference Consortium to have it, or Ensembl, or NCBI, but I could not find it.

Do you mean cytoband coordinates? For the human genome you can find them for example at UCSC.
hgdownload.cse.ucsc.edu/goldenPath/hg38/database/cytoBand.txt.gz

Hey thanks!, do you know why is there in hgdownload.cse.ucsc.edu/goldenPath/hg38/database/ also the file cytoBandIdeo.txt.gz? I checked, and it has exactly the same info as cytoBand.txt.gz. BTW: do you know what the fourth column in cytoBand.txt.gz means? It has values like: gneg, gpos25, gpos50, gpos100, acen, etc..

The answers can be found on google.


text:
Evidence of emergence of "genuine" novel protein interactions by aa mutation
Anyone knows of an example of a protein that, without coming from a recent duplication event, underwent a mutation(s) that caused it to have a novel interaction with a new ligand, substrate, other protein or molecule?
By novel I imply that it doesn't "go back" to an ancestral interaction or doesn’t develops an interaction that it's already performed by a close relative.
I assume that these mutations are extremely rare, and that a duplicate is almost always necessary to have them fixed in a lineage.

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

That's easy, red blood cells and the single mutation causing the switch in charge resulting in sickle cell anaemia. There are five different convergent mutations although the concept of convergence here has been challenged. The mutation results in the physical collapse of the protein structure and massively reduces its binding with oxygen. The receptor interactions with the malarial parasite render infection of the red blood cell impossible. Even in the heterozygote state, oxygen binding with the heme group is not massively disrupted, but the ability of the malaria parasite in particular plasmodium falciparum, to invade the red blood cell is again greatly impaired. It is an incredible example of human evolution and demonstrates the selective pressure of holoedemica malaria. There is no gene duplication involved here and again due stress there are potentially 5 convergent events resulting in the same phenotype and the scale of protein collapse from circular red blood cell to sickle cell shape is the result of what would otherwise appear to be a small change in amino acid charge for a single amino acid substitution.

Thank you, I’m aware of those mutations, but they don’t induce a novel interaction. As you say, they disrupt the tertiary and quaternary structure of hemoglobin. I’m asking for “constructive” mutations. Like for example if some mutation caused hemoglobin to bind preferentially other molecule instead of their known ligands.


text:
NCBI retrieve genome assembly with splicing or "modifications"
I am trying to fetch genome assemblies that include splicing information.
E.g. Here, is a genome assembly of a yeast, where genes show variants like this: join(5928..7626,7685..7785,7857..7982).
I am retrieving the information from NCBI with the following commands:
efetch -db nuccore -id CP003820.1 -format tabular. I tried different formats and it changes but I don't get the join(...)
1) Am I correct and the join(...) indicates "splicing-like" modifications in yeast?
2) How can I get the join(...) (any other format is fine too) with efetch?

The page you link to is a GenBank Flat file, the default format used by GenBank. So, to get that format, which includes the join lines, use:
./efetch -db nuccore -id CP003820.1 -format gb

The join() lines give the different ways in which a sequence can be built from this gene. For example, one mRNA sequence produced by the locus you mention (CP003820) is:
  mRNA            complement(join(100..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..>5645))

The corresponding coding sequence is for the mRNA above is:
 CDS             complement(join(168..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..5645))

These aren't necessarily splicing "modifications", and aren't "variants", they are simply the way that the exons are joined to make an mRNA.
Note that since CP003820.1 is "Cryptococcus neoformans var. grubii H99 chromosome 1, complete sequence", there are multiple genes described in there and each will have its own (or multiple, if there is alternative splicing) join statements. Specifically, I see 1014 genes described in the file:
$ grep -Pc '^\s+gene\s+' CP003820.gb
1014

950 mRNAs:
$ grep -c 'mRNA.*join(' CP003820.gb
950

And 955 proteins:
$ grep -Pc 'protein_id=' CP003820.gb
955


text:
Can we use GTEx data as control data for TCGA data?
I am using Recount2 TCGA data and was wondering is it right to use GTEx data as control data for this. I would really appreciate your views on this?

Definately not directly. Probably not even if used indirectly. Different projects, done at different times, with different protocols, different machines, different library preps etc.
We have tested combining GTEx and TCGA data for several different purposes where one might expect the analyses' internal controls would have allowed compatibilty, but even then we have always found a persistent batch effect which we were unable to remove.

I am not sure how appropriate it is, but GTEx data was used as the normal controls for TCGA data in Comprehensive analysis of normal adjacent to tumor transcriptomes. It was not just a single panel, but a major part of the paper.
Alternatively, there is RNAseqDB:

Here, we developped a pipeline that processes and unifies RNA-seq data
  from different studies. Using the pipeline, we have processed data
  from the GTEx and TCGA and have successfully corrected for
  study-specific biases, allowing comparative analysis across studies.

One thing with RNAseqDB is that they provide ComBat normalized FPKM values, which cannot be used for DGE analysis with methods which require raw counts, such as DESeq2 and similar.


text:
How to merge .fastq.qz files into a single .fastq.gz with their same id without losing any content in parallel
I have a large number of .fastq.gz files of different lane and reads. I have to merge them each reads group files into single .fastq.gz files.
**eg:
1st type
NA24694_GCCAAT_L001_R1_001.fastq.gz
NA24694_GCCAAT_L001_R1_002.fastq.gz
…....
…....
NA24694_GCCAAT_L001_R1_033.fastq.gz
2nd type
NA24694_GCCAAT_L001_R2_001.fastq.gz
NA24694_GCCAAT_L001_R2_002.fastq.gz
…....
…....
NA24694_GCCAAT_L001_R2_027.fastq.gz
3rd type
NA24694_GCCAAT_L002_R1_001.fastq.g
NA24694_GCCAAT_L002_R1_002.fastq.g
…....
…....
NA24694_GCCAAT_L002_R1_040.fastq.g
4th type
NA24694_GCCAAT_L002_R2_001.fastq.gz
NA24694_GCCAAT_L002_R2_002.fastq.gz
…....
…....
NA24694_GCCAAT_L002_R2_040.fastq.gz
so now i need to merge all the corresponding files into single file.
Output:
EA00694_GCCAAT_L001_R1.fastq.gz
EA00694_GCCAAT_L001_R2.fastq.gz
EA00694_GCCAAT_L002_R1.fastq.gz
EA00694_GCCAAT_L002_R2.fastq.gz
I tried ""cat'" but its not possible to do for large files. Could any one help me to run in script parallel in perl or python or shell. or any linux command which can be used for large file.
Thanks all.

Why is cat not possible for large files, please tell more about the size and number of files and/or the error you get when running the cat command.

the folder size 40gb.. there are >200 files for each type.. so I couldn't do it manually thats why I am looking for scripts or any efficient command

Did you try cat NA24694_GCCAAT_L001_R1_???.fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz?

You need to dump your data on an external harddrive with sufficient. If you run out of harddrive it will throw an exception, or destroy the data. You are requesting a solution in higher order programming, I don't think its needed, just carefully crafted key words to build the file up in stages rather than all at once.

All you need is cat. You won't find any better tool for a simple job like this. Just run:
cat NA24694_GCCAAT_L001_R1*fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz
cat NA24694_GCCAAT_L001_R2*fastq.gz > EA00694_GCCAAT_L001_R2.fastq.gz
cat NA24694_GCCAAT_L001_R3*fastq.gz > EA00694_GCCAAT_L002_R1.fastq.gz
cat NA24694_GCCAAT_L001_R4*fastq.gz > EA00694_GCCAAT_L002_R2.fastq.gz

That should work just fine. You say you tried cat and it "didn't work" but since you don't tell us how it failed, I can't really help. The only issues I can think of is that either there are too many files (as in several hundred thousand, whatever the value returned by getconf ARG_MAX on your system which is 2097152 on mine) or, more likely, you are running out of disk space.
If it's a disk space issue, you might be able to get around it by adding each file and then deleting it:
for i in 1 2 3 4; do
    for file in NA24694_GCCAAT_L001_R${i}_*fastq.gz; do
        cat "$file" >> EA00694_GCCAAT_L001_R${i}.fastq.gz && rm "$file"
    done
done

IMPORTANT: the command above is destructive. It will delete each file after it has been added to the new one. If everything works fine that's not a problem, but if I've made a mistake or your file names are slightly different you might lose your data. So I strongly urge you to simply run the cat commands in the first section on a machine with enough disk space to deal with them.

I was using cat command manually by adding to file names to it  I had these types of files    EA13698_AGTAAT_L001_R1_001.fastq.gz
EA13698_AGTAAT_L001_R2_032.fastq.gz
EA13698_AGTAAT_L002_R1_003.fastq.gz
EA13698_AGTAAT_L002_R2_037.fastq.gz
I had these the four types of files, in each group it had more than 50 files so I couldnt keep on adding names to cat command.   eg "cat EA13698_AGTAAT_L001_R1_001.fastq.gz
EA13698_AGTAAT_L001_R1_002.fastq.gz
.... > EA13698_AGTAAT_L001_R1.fastq.gz. manually doing this it takes time so I was searching for script which will do looping.

No need to loop, just do them all at once with globbing. 50 files won't be too long and you won't have to type them out.

IMO the way I would do it is with ">>" (append) rather than brute forcing the concatenation with ">" (over-ride) ... which is what the loop is about. I like the shell loop.

As said by @terdon, I had used * to group the types of file in 'cat'.. it worked. Thanks all.

We sequence DNA on daily basis. Like terdon said, you can use cat to group all the FastQ of all your sample whether they are bgzipped or not like so:
root@slurm1: ~  # ls
am224_S6_L001_R1_001.fastq.gz
am224_S6_L001_R2_002.fastq.gz
am224_S6_L002_R1_003.fastq.gz
am224_S6_L002_R2_004.fastq.gz
am232_S10_L001_R1_001.fastq.gz
am232_S10_L001_R2_002.fastq.gz
am232_S10_L002_R1_003.fastq.gz
am232_S10_L002_R2_004.fastq.gz

You could do:
# initiate am empty array
listSamples=()

# loop to add the unique samples to an array
for i in *.fastq.gz;
do
    if [[ "$sampleName" =~ "${listSamples[@]}" ]]; then
        sampleName=$(echo $i | cut -d'_' -f1)
        listSamples+=("$sampleName")
    fi
done

# loop in the newly created array of samples and group them by reads number
for i in "${listSamples[@]}"
do
    cat "$i"_*_R1_* > "$i"_1.fastq.gz
    cat "$i"_*_R2_* > "$i"_2.fastq.gz
done

However, I highly suggest you to compress your files as .bz2 to save storage as much as possible. You can achieve this using the multithreaded tool lbzip2. Assuming that you have 10 CPUs available, you can do:
for i in "${listSamples[@]}"
do
    bzcat "$i"_*_R1_* | lbzip2 -n 10 --best > "$i"_1.fastq.bz2
    bzcat "$i"_*_R2_* | lbzip2 -n 10 --best > "$i"_2.fastq.bz2
done

Hope this helps !

I like the bz2 storage suggestion


text:
How can I calculate coverage at single bases using a bam file?
I'm looking for a way to input a vcf or bed file (with specific base positions) and a bam file, and get the coverage at each base position (ie single base bins) using the bam file. I also want the strand information so ideal output would be something like:
chr base_position total_coverage fwd_coverage rev_coverage

I've tried bamCoverage and bedtools Coverage but with bamCoverage, I can't figure out how to get single-base pair resolution and with bedtools Coverage, I can't figure out how to get the strand coverage information.

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.

git clone https://github.com/lh3/htsbox
cd htsbox && make
./htsbox pileup -cCf ref.fa aln.bam | less -S

This output a VCF containing positions covered by at least one read. Tag ADF and ADR give the forward/reverse depth for each allele. You can combine them to get the total counts on both strands.

Brent Pedersen (@brentp) claims mosdepth is 2x as fast as samtools depth. I'm not sure whether this provides strand information, but it might be the fastest tool available.
Another tool of possible interest (also written by Brent) is indexcov. This doesn't give single bp precision, but is able to provide even quicker estimates of coverage at 16Kb resolution.
Installing these tools from source is non-trivial, and I struggled with installation in the past. However, they are both now available in bioconda, which should make installation simple.


text:
Tabix queries spanning chromosomes
I can do something like this:
$ tabix ${url} ${range}

where ${range} spans a region on a single chromosome, e.g., chrA:X-Y.
Can a tabix query span multiple chromosomes, e.g., where ${range} looks like chrA:X-chrN:Y (or similar), which returns data between chrA:X and chrN:Y (using either a lexicographical or natural sort order between chrA and chrN)?
Note: I do not want to pass in a BED file of intervals. If it can't be done without passing in a separate file, that's fine. I only want to provide a one-liner to give to someone who can run tabix on their end, to do this type of query directly and remotely.

Looking at the htslib source code, it seems like the answer is no, a single region specifier cannot span multliple chromosomes.
Seems like the closest you can get is to enumerate all the chromosomes as separate regions on the command line. So to get everything between chrA:X and chrN:Y, you would have to write:
tabix inputfile.gz chrA:X chrB chrC ... chrN:1-Y

(If you omit the end coordinate like chrA:X, the region encompasses chrA from the specified coordinate to the end. If only specify the chromosome name without the coordinates, like chrB, the region encompasses the whole named chromosome, and for the last chromosome, you can specify up to the desired point by starting from the beginning: chrN:1-Y.)

Okay, thanks, that worked great. I was able to order chromosomes as desired to write standard output in BEDOPS sort-bed order.


text:
Queries regarding MetagenomeSeq code
I am trying to understand metagenomeSeq, I will be very much thankful if you could help me to understand this:
I do have three treatment condition (T1, T2, and T3) and want to see how these treatments affect the microbial diversity in rice roots. I am interested in the pairwise comparison between them.
I used below code but I am not sure whether it is right or wrong, I do also have few queries regarding below code. I initially imported whole data into phyloseq and then filter the root data and convert into metagenomeSeq object.
root_data <- phyloseq_to_metagenomeSeq(root_M_filter)

DO I need to used normalize count for further processing? what is setting?
normalizedMatrix <- MRcounts(root_data, norm = TRUE)
settings = zigControl(maxit = 1, verbose = FALSE)

Is it correct way to make model matrix?
mod = model.matrix(~root_data$Treatment)
colnames(mod) = levels(root_data$Treatment)
res = fitZig(obj = root_data, mod = mod, control = settings)
res
An object of class "fitZigResults"
Slot "fit":
An object of class "MArrayLM"
$coefficients
                                      T1           T2          T3
71f84e7910006f22684121564206e8ca -2.8129112 -0.321366092 0.06113294
54e204fb99c80764e964456dadd6a0e5 0.5643061 0.019018800 1.01351768
55cd1fc570879d645bbf7a3642e9b0a8 -3.1159264 0.072348076 0.02205983
65f5c31e12c277aec319e2096463f9d2 -0.1355134 -0.004567206 0.53119153
7bed62f0fef250fd831dcf13bf43f4fc 1.9469670 -0.596315540 -0.35124338
                             scalingFactor
71f84e7910006f22684121564206e8ca 1.62101725
54e204fb99c80764e964456dadd6a0e5 -0.24292482
55cd1fc570879d645bbf7a3642e9b0a8 1.50462644
65f5c31e12c277aec319e2096463f9d2 0.05833637
7bed62f0fef250fd831dcf13bf43f4fc -0.56247683

What is design here?
finalMod = slot(res, "fit")$design

It is correct option for contract.make, or should I also involved some
more options?
contrast.matrix = makeContrasts(T1 - T2, levels = finalMod)
fit2 = contrasts.fit(zigFit, contrast.matrix)
fit2 = eBayes(fit2)
topTable(fit2)
                                 logFC   AveExpr         t      P.Value
3147790f0d5a78316fb9dd64f53b9473 11.262064 14.091887 16.331949
3.482815e-15
d2ec9f3b77975c0f457e4b7413b217ff 9.466087 12.212476 10.606122
6.159311e-11
7b908379d1fd2e75d84fd9aaecfd3f77 6.907182 9.784971 9.163780
1.267984e-09
d5fdc1df6dfd31acce9a34b9d8a76590 5.404459 8.373729 7.035926
1.802316e-07
19ddf8e59673c884734cf6093c672da0 5.381863 8.449413 7.011910
1.912363e-07

You get warnings, not errors. The problem is you keep trying to load ph_mgs as a dataset from a package, while it is an object you create. However, this warnings do not prevent you to analyze the ph_mgs as you wish

@DrYogesh If the problem is unrelated then I would recommend deleting this post. But you could also [edit] this one to include the relevant section (if your problem is highly related to this code).

Please use code formatting not quotation formatting for the question.


text:
Alternative splicing positions on protein
I have mass-spectrometry data (Human cellsanalyzed and annotated by MaxQuant).
I want to find which peptides are covering positions of former splice junctions (on the mRNA sequence), based on protein id and peptide positions.
Where could I find an appropriate data? I want to know, based on peptide-sequence, annotaion and position on protein, whether this peptide is a single exon, or if it contained areas of former splice junctions.
Itai

Hello Itai, what do you mean by "appropriate data'? Do you mean where you can find annotations of splice junctions? What species? Please, [edit] your question and be clear and verbose. The second thing that could be improved is that your question shows very little research effort. What have you tried? Have you found any clue how the problem could be resolved? Yes? Write it here.

I edited the question

Well done, It's way clearer now. Maybe just add the species you are looking at and the results of your own research on the topic...

If you are considering the human genome as "data" (not mRNA), I would run your single peptide via tblastn which will take your protein(s) translate them to nucleotides and run a blast against Genbank. The splice sites will be immediately apparent because they force indels in your query. If there are no indels there are no introns (splice sites). If you want to process several proteins, one solution is to use eblast, a command line version of tblastn  (Anaconda package). You can confirm this against a full genome annotation (note, I don't work on humans).
Alot of work has been done on splice sites, a friend worked on this, but I think you just want to focus on a restricted number of target peptides in your mass spec. output.
mRNA should be avoided IMO if you are looking for splice sites, i.e. intron boundaries.

eblast sounds like a great idea, thanks!