text:
awk for motif splicing from an amino acid alignment
I'm looking to splice out specific amino acid alignment positions as a continuous motif from an alignment (~15 000 sequences). The code is written in Perl (I agree not many developers now) but splicing method is awk which a lot of you know. Please note the hash is "tied" (not shown), this just means the order of the sequences is preserved.  I can get the awk code below to return single amino acid positions.
I am struggling to get a loop working in awk.
I agree the Perl wrapper makes it more complicated, but its part of a huge data pipeline (all written in Perl).
#!/usr/local/ActivePerl-5.20/bin/env perl
#use strict;
use warnings;

# I loaded the alignment into a hash, where key = sequence and value = protein
my %dna = (
    virus1 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF",
    virus2 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF",
    virus3 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF"
);

# I define the start and end of the alignment positions I am splicing
my $n = 1;
my $z = 5;

# I turn the alignment into a string with separated by carriage returns
my @array = map {"$dna{$_}"} keys %dna;
my $string = join "\n", @array;

# This is the awk loop, note in Perl the "``" means issue a shell command

my $splice = `echo "$string" | awk 'BEGIN {FS=""} {x = $n; while (x < $z); print \$ x; ++x}' 2>1`;
print $splice;
__DATA__

Expected output
LGVLV
LGVLV
LGVLV

Observed output
I get the original sequences returned
I am familiar with Python/Perl solutions of line by line processing but awk is more efficient IMO.

Could you [edit] and show the output you're expecting? I speak both Perl and awk, but I don't get what output you want. Also, I really doubt that calling awk will be more efficient than doing it natively in perl. You don't seem to be using anything that awk is particularly good at.

Amended, thanks. I will post the solution I got working, but this is because I don't understand awk

The main problem is that your while only encloses the first statement after it, and not the second. You need { } for that. If you could explain what that awk is supposed to be doing (I'm pretty sure you've copied it wrong, is there really a space between \$ and x?), I should be able to post an answer. Do you just want the first 5 characters of each sequence? Don't you want positions $n to $z (3-7)?

I think you are trying to use awk to get characters 3 through 7 (inclusive) of each input sequence. If so, a working version of your approach would be this:
my $splice = `echo "$string" | awk 'BEGIN {FS=\"\"} {x = $n; while (x < $z){ printf \"%s\",\$x; ++x} print ""}' 2>1`;

That, when added to your perl script, produces:
$ foo.pl
VLVIL
VLVIL
VLVIL

That isn't the output you show, but the output you show is just the first 5 letters of each sequence, which doesn't seem to correspond to what your script is doing.
Now, the reason your version didn't work, is that the while was like so:
while (x < $z); print \$ x; ++x}

Assuming the \$ x was a typo, the problem is that the ++x isn't part of the while loop. So that should never even have exited. Seeing as we're in the same time zone, and given the time, I am guessing you're tired and that you miscopied something when trying to give us a reproducible example. In any case, to get the loop right, you need {}:
while (x < $z){print \$x; ++x}

However, if all you want is a substring, a loop is a very inefficient way of doing that. You can instead use awk's substr command:
  substr(s, m[, n ])
             Return  the at most n-character substring of s that begins at position
             m, numbering from 1. If n is omitted, or if n specifies  more  charac‐
             ters than are left in the string, the length of the substring shall be
             limited by the length of the string s.

So, in your script, that would be:
my $splice = `echo "$string" | awk '{print substr(\$0,$n,5)}'

Since you're already using a language that excels in string manipulation, I really don't understand why you would try to use awk instead. It will not be faster or more efficient than the same thing in perl:
print substr($_,$n,5) . "\n" for @array;

By calling awk you are also calling the shell which will be slower by definition, especially if you're echoing thousands of lines. In fact, I tried this with 15000 sequences (just repetitions of your first one) and using either the awk loop I gave above, or the awk+substr, or the native perl substr:

awk loop, average time taken to run the same thing 100 times: 0.0525
awk substr, average time taken to run the same thing 100 times: 0.0522s
perl substr, average time taken to run the same thing 100 times: 0.0464

And you can make it even faster by skipping the array and iterating over the keys directly:
print substr($dna{$_},$n,5) . "\n" for keys(%dna);

That averaged at 0.0376s over 100 runs.

Huge thanks @terdon. You are right I should have performed a runtime check :-( such as use Timer::Runtime. I worked on the implicit assumption awk was faster for this operation. David Cross "lectured" on runtime in detail in "Data Munging for Perl": There is no denying that 0.0525 vs. 0.0376s is 30% faster.

@MichaelG. in my experience, awk is blindingly fast when your data are organized in fields because of its implicit splitting. If not, however, it doesn't give any real advantage over a language like Perl which is simply brilliant for text parsing.

If you have the sequences in a seperate file like this:
$ cat input
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF

you can use a simple cut to fetch a substring:
$ cut -c1-5 input
LGVLV
LGVLV
LGVLV


text:
Multi-line efetch digging
I'm using efetch to dig out a middle "envelope" protein within a viral genome specifically within fasta format. The code below works but may not be robust for an entire database. Thus I want to replace the two "perl" lines with formal code within the Entrez package. I tried using efilter -query but couldn't get it to work. Genome positions are one approach but are not absolute from one virus to another. Rewriting this may be the best option.
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" |\
xtract -pattern INSDSeq -element INSDQualifier_value |\
perl -p -e 's/^(.+?)\s+taxon.*.*\s+\b([A-Z]{300}[A-Z]+).*/>$1\n$2/'|\
perl -p -e "s/[\t \;\/]+/_/g"

Output
>InfluenzaAvirus_A_PuertoRico_8_1934
MDVNPTLLFLKVPAQNAISTTFPYTGD.....

I will provide a more detailed explanation shortly.
Any ideas welcome.

In the meantime, you could streamline the perl a little to | perl -pe 's/^(.+?)\s+taxon.*envelope protein.*\s+\b([A-Z]{300}[A-Z]+).*nonstructural protein 1.*/>$1\n$2/; s/[\t ;\/]+/_/g. There's no need to process the file twice and you also don't need to escape the ; inside a character class.

Hi @terdon, yes and thanks now you've mentioned it, it's obvious.

Not exactly what you were looking for but may be you can use something like this? efetch -db protein -format gpc -id "$ncbi2" | xtract -pattern INSDSeq_feature-table -group INSDFeature -if INSDFeature_key -equals 'mat_peptide' -tab '\n' -element INSDFeature_location,INSDQualifier_value | grep 'envelope'

Perfect it works a treat @vkkodali. Huge thanks. A small amount of regex will finish the fasta format. Maybe you would like to provide a formal answer so I can sign this question as complete?

Not exactly what you were looking for but may be you can use something like this?
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" \
    | xtract -pattern INSDSeq_feature-table \
        -group INSDFeature \
        -if INSDFeature_key -equals 'mat_peptide' \
        -tab '\n' \
        -element INSDFeature_location,INSDQualifier_value \
    | grep 'Segment 2'

For explorations of the INSDSet XMLs, xtract has -insd option that you can use. In this particular case, you can do something like this:
efetch -db nuccore -id "$ncbi2" -format gpc \
    | xtract -insd mat_peptide product peptide

For more information, checkout EDirect documentation and some  example queries

Thank you its perfect. Stacked together in a small subroutine and parsed with a bit of regex its fantastic. "esearch -query" is used to obtain the accession numbers.

source: Stack Exchange

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

sort: votes


text:
PDB file downloading: pymol automation vs. manual
I automated a PDB download using a Pymol script (below)
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x)
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

When the script is run within Pymol it pulls down a collection of pdbs, however this flopped when I  plugged it through a Biopython PDB parser, but was okay when I manually downloaded the same pdb file from RCSB.org
A manual downloaded PDB was perfect, meaning the amino acids which cannot be observed by X-ray crystalography (e.g. they not sufficiently stationary) were recognised as being absent from the parser (resulting in a protein ~390 amino acids), but the Pymol scripted download contained all amino acid residues, i.e. even those which couldn't be observed in the X-ray (resulting in a protein of ~440 amino acid residues). Any ideas why the script resulted in a pdb that flopped and how to correct it?

In PyMol 1.8+ fetch downloads mmcif files by default. Your script must convert mmCIF to pdb, that's why the result differs from a pdb file downloaded directly from wwPDB.

To expand on the comment by marcin:

fetch downloads files in mmCIF format by default (https://pymolwiki.org/index.php/Fetch). Not all PDB entries have PDB format files, e.g. due to too many chains. Presumably this is why the change was made, though mmCIF files tend to be larger and hence download slower.
When calling save with the .pdb file extension the structure is converted to PDB format. Perhaps the converter uses the SEQRES records to define the sequence, and hence gives bad files when there are some residues not resolved.

To correct it you could try:

Saving as .cif and reading in with the Biopython mmCIF parser.
Providing the type=pdb argument to fetch, which will force PDB format downloads. Some PDB files may not be available.
Don't do it in PyMol and instead use the Biopython structure downloader (http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc187).

Also, if download speed is limiting you then consider using the MMTF file format, which is a binary format and hence smaller. Biopython can read in and write out MMTF files.

Thanks point 2 was exactly how I solved the issue (and it worked). Let me look at the Biopython solution. I looked in detail at MMTF but it that appeared more aimed at pdb->MMTF.

So this was the code I used and it worked (as far as I could see) and refers to point 2 of @jgreener's advice and was based on the comment by @marcin
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x, type='pdb')
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

There are some Bio SE members who have serious skills in Pymol scripting, and constructive critique would be very welcome.

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


text:
efetch nucleotide -> protein ids
I'm downloading translated viral genomes initially via Blast - which had shortcomings - and now by efetch.
What I want to do is obtain a complete list of all translated protein ids from a nucleotide Genbank entry.
efetch -db nucleotide -format gpc -id MT039887 | xtract -insd ....

The protein id and the translated gene sequences are within the efetch -db nucleotide -format gpc -id MT039887 output ... but how do I use xtract to parse them?

If you want the cds in fasta format, try:
esearch -db nucleotide -query MT039887 | efetch -format fasta_cds_na

And for the proteins sequence in fasta format:
elink -db nucleotide -target protein -id MT039887 | efetch -format fasta

Wow!! and thank you. Thats the answer ... days like this BioSE works. @zorbax I owe you. "Acknowledgements" thank you heading your way (thats not a metaphor).

In case you are still interested in just getting the protein ids and not sequence, you can use xtract like this:
efetch -db nucleotide -format gpc -id MT039887 | xtract  -insd CDS gene product protein_id
MT039887.1  orf1ab  orf1ab polyprotein           QHZ00388.1
MT039887.1  S       surface glycoprotein         QHZ00389.1
MT039887.1  orf3a   orf3a protein                QHZ00390.1
MT039887.1  E       envelope protein             QHZ00391.1
MT039887.1  M       membrane glycoprotein        QHZ00392.1
MT039887.1  orf6    orf6 protein                 QHZ00393.1
MT039887.1  orf7a   orf7a protein                QHZ00394.1
MT039887.1  orf8    orf8 protein                 QHZ00395.1
MT039887.1  N       nucleocapsid phosphoprotein  QHZ00396.1
MT039887.1  orf10   orf10 protein                QHZ00397.1

Thank you, this is really useful.


text:
ETE4 manual or improved code description
ETE3 * (Environment for Tree Exploration) is a major package and in particular Python package to navigate phylogenetic trees. It is the premier coding library/package for tree manipulation.
ETE4.0.0 beta has been released via GitHub. The download contains an old ETE manual (2009) and a README file for ETE4, thats about it. Does anyone know the location of the new ETE4 manual, or list of functionality beyond what is in the README (posted below). Reading the code is the last resort ;-)

README ETE4

Read and write support for trees in Newick format
Multiple functions for traversing, searching and manipulating tree topologies and node annotations
Integration with NCBI Taxonomic database
Integration with GTDB database
Programmatic visualization framework
Visualization of huge trees based on adaptive zooming
Comparing trees
Phylogenetic functions:

orthology detection
phylogenetic distance?

Reference
* ETE 3: Reconstruction, analysis and visualization of phylogenomic data.
Jaime Huerta-Cepas, Francois Serra and Peer Bork.
Mol Biol Evol 2016; doi: 10.1093/molbev/msw046

The github page you link to says they have documentation on their official website (http://etetoolkit.org/) and, indeed, clicking on the documentation link there took me to  http://etetoolkit.org/docs/latest/index.html. Is that what you're looking for?

Thanks @terdon, this is the documentation for ETE Toolkit 3.0.0b34 (last updated Mar 18, 2016). I recognise the ETE3 to ETE4 is small outside phylogenetics, but it is a big upgrade. The authors may respond here, if not I'll send them this link.

I suspect you will get more traction if you just ask on their github. This isn't a question that is likely to be answerable by anyone other than the devs.

Good point. https://github.com/etetoolkit/ete/issues/651 ... will pop their response as an answer here. ETE4 has potential.

Just to feedback, the GitHub ticket was:

ETE4.0.0 beta contains an old ETE manual (2009) and a README file for ETE4. Is there an updated ETE4 manual beyond etetoolkit.org/docs/latest/index.html or the README of ETE4?

lm-jkominek wrote,

I know it's a bit of necro, but I wanted to second this request  - while the current docs work fine enough for ETE3, they have diverged quite far from ETE4 and there's a bunch of errors that come through when using the latest version that have no explanation...

Thus I'm not alone ... (nothing to do with UAPs BTW)
jhcepas wrote,

absolutely. We are finishing some final changes to the core of ETE4 and will start with the documentation. In the mean time, @jordibc has a list of most important changes regarding compatibility. Could you share somewhere as a wiki page @jordibc ?

jordibc replied,

Yes, we are working on an updated version of the manual to be released with ete4. In the meantime, most of the old code can be run with some of the changes mentioned in the part of the wiki on migrating from ete3 to ete4. @lm-jkominek and @sarscov2 thanks and hope this helps for the moment!

The GitHub wiki was updated here. The Wiki page is kept up to date,

"Jordi Burguet-Castell edited this page 4 days ago · 2 revisions"

The actual additional features (why ETE4 is needed) and how to access them are not listed here. However the changes in syntax are listed. Examples of changes ...
# ETE3
t = Tree(filename)
# ETE4
t = Tree(open(filename))
# ETE3
Tree(name='Hi', dist=1, support=1)
# ETE4
Tree({'name': 'Hi', 'dist': 1, 'support': 1})
# ETE3
Tree(nw, format=1)
# ETE4
Tree(nw, parser=1)

In summary there's different syntax. So we can start using ETE4

thank you very much M__


text:
Sequence format: grouping sequences within a flat file
Is there a sequence format allowing multiple sequences to be grouped  across a flat file? The specific application is stacking separate nucleotide sequences against a given RNA secondary structure.
To rephrase the same question I'm looking to stack viral nucleotide data in different groups across a single flat file. FastaQ provides an example:
Traditional fastQ,
@sequence_name
AGCTGAT
+
...(quality score)

@sequence2
AGCTGAC
+
...(quality score)

What I'm looking for is a format,
@Gp1_virus1-2-3
AGACGTAG
+
AGAGACAA
+
AGAGACAT

@Gp2_virus1-2
AGATGTAG
+
AGATACAA

In the above example there are two groups comprising 3 (group 1) and 2 (group 2)  individual sequences (viruses). Here 'group' represents different structures.
Stockholme format has potential. My understanding of this format, which is basic, is it stacks a sequence, together with structural information. What is sought is stacking multiple sequences and then a structure format. Any clarification would help.
There are multiple XML formats, which may be relevant.
Just to mention a 'flat file' is needed to allow manual inspection. Dumping JSON is a solution, but isn't friendly for optical inspection.

Application
The immediate downstream application is structure/sequence viewers. I think at least one. viewer will accept Stockholme 1.0 format. rFAM is very similar (@terdon's comment)
Summarised Question The question could be rephrased as to whether multiple RNA families can be stacked inside one Stockholme 1.0 or else rFAM format?
Some of the members here are likely familiar with pFAM. So could multiple protein structures be stored and retrieved inside one pFAM file.
@gringer's contribution is very welcome.

What is this for? I mean, since you don't seem to need something that will be understood by a specific downstream program, can't you just define your own format?

Thanks, @terdon I've updated the post. A hybrid solution may be needed (last resort).

An unmapped SAM file format allows for read groups. You can use samtools import to create unmapped SAM files from FASTQ files:
samtools import -R virus-1-2 in12.fq > 12_reads.sam
samtools import -R virus-1-2-3 in123.fq > 123_reads.sam
samtools merge out.sam 12_reads.sam 123_reads.sam

Perfect, thank you. It can be manipulated via Biopython SamtoolsCalmdCommandline

I'm awarding 'solved' to @gringer, because thats the context of the question.
The best solution I could find is stacking Stockholm (sth)  formats together. In this example the RNA families are associated with one structure:
# STOCKHOLM 1.0
#=GF ID    UPSK
AF035635.1/619-641             UGAGUUCUCGAUCUCUAAAAUCG
M24804.1/82-104                UGAGUUCUCUAUCUCUAAAAUCG
J04373.1/6212-6234             UAAGUUCUCGAUCUUUAAAAUCG
M24803.1/1-23                  UAAGUUCUCGAUCUCUAAAAUCG
#=GC SS_cons                   .AAA....<<<<aaa....>>>>
//

Biopython allows sth format to be input.
from Bio import AlignIO
from pathlib import Path

file = 'test.sth'
path = '/Volumes/mypath/testSth'
alignment = AlignIO.read(open(Path(path, file)), "stockholm")
print("Alignment length %i" % alignment.get_alignment_length())
print (alignment.column_annotations['secondary_structure'])

Output,
Alignment length 23
.AAA....<<<<aaa....>>>> # this is the RNA structure

If we look the alignment object via
from pprint import pprint
pprint(vars(alignment)) # or via print(alignment.__dict__)

The output

{'_annotations': {},  '_per_col_annotations': {'secondary_structure':
'.AAA....<<<<aaa....>>>>'},
'_records':
[SeqRecord(seq=Seq('UGAGUUCUCGAUCUCUAAAAUCG'),
id='AF035635.1/619-641', name='AF035635.1',
description='AF035635.1/619-641', dbxrefs=[]),
SeqRecord(seq=Seq('UGAGUUCUCUAUCUCUAAAAUCG'), id='M24804.1/82-104', name='M24804.1', description='M24804.1/82-104',
dbxrefs=[]),
SeqRecord(seq=Seq('UAAGUUCUCGAUCUUUAAAAUCG'), id='J04373.1/6212-6234', name='J04373.1',
description='J04373.1/6212-6234', dbxrefs=[]),
SeqRecord(seq=Seq('UAAGUUCUCGAUCUCUAAAAUCG'), id='M24803.1/1-23', name='M24803.1', description='M24803.1/1-23',
dbxrefs=[])],  'annotations': {}} {'secondary_structure':
'.AAA....<<<<aaa....>>>>'}

The reason vars(alignment) is required is to internally construct an alignment object in Stockholm format, here's the first part of the skeleton 'pseudo-code',
from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord

outAlign = MultipleSeqAlignment([])
record in alignment:
    name = record.id
    myseq = record.seq
    sr = SeqRecord(Seq(name), myseq, '', '')
    outAlign.append(sr)

The above code constructs a basic alignment, the bit thats missing is to add the following key/value to the alignment object of {'secondary_structure': '.AAA....<<<<aaa....>>>>'}.
Its complicated, but you get use to it.

https://github.com/biopython/biopython/blob/master/Bio/AlignIO/StockholmIO.py


text:
Data sets for prostate/cervical cancer
I am seeking data sets for prostate/cervical cancer. I'm investigating how close synthetic data (genAI) can approximate field data.
In the case of prostate cancer I'm looking for:

Age
PSA
Gleason
Stage
Treatment
Side Effects
Outcome

65
10.5
7
II
Radiation
Fatigue
Stable

I have investigated several repos for the data and numerous papers for example. The National Cancer Institute (first hyperlink) requires a formal request and these can take a long time and it's not clear whether the resulting data is appropriate. The data set of the paper isn't available but the code is available (I am not sure Elsevier demand it, Nature and PLoS journals do). In theory I can write to the authors and then the editor if that fails. Again this takes a long time but it depends on journal rules
Again, if anyone does have experience with this data set type, it would be really helpful, because even if I have to submit a formal request I can be confident the data will be what is sought.

Just to report back, I found a reasonable data set on Kaggle:
https://www.kaggle.com/datasets/ukveteran/stage-c-prostate-cancer

The top two lines of the data set are posted above, where PG is Proliferation Growth,
G2 is cell cycle G2 phase, ploidy - we all know, EET is  Epoxyeicosatrienoic Acids which is a topical biomarker. It is not a huge data set at 146 patients, but it's enough. I have a fair amount of expertise in genAI. Note, this is an image because the original source data is under Kaggle ownership.
In case you are wondering generative AI (gemAI) is in its infancy in the biotech sector, but that is likely to change and the ability of the largest model to have ingested key data sets in their training is of huge interest. That can only be determined empirically. To give some idea, protein generative AI was part of the last Nobel prize for David Baker. Thats how disruptive its going to become.
There are a fair few miRNA studies on prostate cancer, so there's a major research front in this area.
A second dataset is here:
https://www.kaggle.com/datasets/nikhilvishwanath/prostatecancer

Logarithmic Prostate Weight is weight
LCAVOL typically refers to "Logarithmic Cancer Volume
LBPH often refers to "Lobular Benign Prostatic Hyperplasia
LPSA typically refers to Logarithmic Prostate-Specific Antigen
Locally Confined Prostate cancer (LCP)
Proliferative Growth possibly pgg45 (dunno)
SVI in the context of prostate cancer typically stands for "Seminal Vesicle Invasion

I cross posted on Biostars later because its a time sensitive request here:
https://www.biostars.org/p/9605437/
There was zero response.


text:
API/website for blasting peptide sequence against database of all plant proteomes
I would like to blast a peptide sequence against all available plant proteome databases, i.e. a blastp metaproteome analysis.
What tools are available? Which database can I use?

Welcome to the site lars. You seem aware of blastp, did it not solve your problem? Or do you expect something else of blastp? Could you explain what have you tried?  (to avoid suggesting the same you have already done)

Thank you for the very quick reply. So far I installed ncbi-blast-2.7.1+. With makeblastdb I generated *.phr, *.pin, *psq from a fasta. I now try blastp against this database. Next steps would be to compile a metadatabase from all plant proteome databases I can find.
I wonder if all of this work has already been done. Is there a website/API that I can use? I simply want to find all plant proteins which contain the amino acid sequence TEKKRRETVEREKE. Thanks.

Please [edit] the question to include this information. In addition, what sequences are in your fasta file? I think that you can download fasta files from NCBI. You  could download the sequences from plants proteins and make the database from that (if this is not already what you did)

At the moment, I work only with arabidopsis thaliana. Is there a plant meta-database? My main objective is to figure out in which plant species this peptide is expressed.

For example, PPDB is only for arabidopsis thaliana and maize.

Sorry I don't work with plants (nor with blast) and I can't answer where to find such data. But add this relevant information to the question, so make it easier to those who might be able to answer you.

A simple BLAST on the UniProt website does the trick (target database plants).
Best matches are subsequences ETVEREKE and EKKRRETV in Green alga and Eelgrass respectively.


text:
How to get PSSM of protein by protein sequence with python?
I didn't find way that convert protein sequences in  FASTA to PSSM by python,
and how to write program in python to find PSSM for our data set ?
Maybe better to question how to find PSSM from NCBI alignments ?
Thank you

Welcome to the site. Where did you search? Did you look up to biopython? What is your data set? A collection of fasta sequences or alignments?

Thank you for your answer , yes i check biopython , but there is no method or class that get protein sequences in fasta and return pssm  . for example from NCBI site , for each sequences.

Ok, so you want to retrieve from the NCBI site, download the protein sequences and then create a pssm? That is different from what you posted. Please [edit] the question to include this details. Would you be fine with a solution that worked with a local blast program?

I think better to use psi-blast from NCBI to find my alignments for each sequence  i don't know easy way to  get aligments and pass them to function to get PSSM (i dont know it is necessary to find pssm)

Sorry I don't understand your last comment. What are you trying to do/get? To get a PSSM you need an alignment of sequences. So you need several sequences and then align them (that's what psi-blast does). But you seem to ask about retrieving the PSSM matrix from a psi-blast search. Have you tried something else like hmmer?

Sorry , i am too beginner . Can you give me  a solution to get PSSM for each of my protein sequence with use psi blast NCBI ? I can calculate PSSM localy on my computer but how to get aligmnets?

Sorry it is not easy to answer your question. The alignments must be done downloading all the sequences of interest and aligning them via t-coffee or similar tools. Then you can calculate the PSSM matrix (not the other way round)

@Llopis Thanks, do you know how to automatically download alignments for each protein sequence ?

Download from where? BLAST perhaps? In what format, and in what language ?

From NCBI , PSI-BLAST, in xml fomat, in python

This would require a major re[edit] of the question. Could you please update the question to keep the question and answer related?

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


text:
Double-counting coverage of overlapped read pairs
EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region.
I used samtools depth to calculate the depth of coverage for samples in the whole Exome region using a GRCh37_ref.bed. These samples are sorted and duplicate marked. I ran this calculation on a few hundred samples to determine how much more sequencing needed to be done and found something interesting. Some samples had 100X coverage at 96,000,000 mapped reads and some samples had 100X coverage at 83,000,000 mapped reads. The samples that had 100X coverage at 83,000,000 reads had read pairs overlapping certain regions of the bedfile (read 1 was covering the same coordinates as read2 to some extent).
I then found out that samtools depth double counts these overlapped regions even though they are technically from the same molecule in sequencing and would be a duplicate bases in the read. See the image for an example:

Each base in this Duplicate Overlap region is double counted when it shouldn't be. How do you guys deal with this duplicate overlap region when calculating coverage across a bam file?
samtools (v1.7):
samtools depth -a -b exome_targets_GRCh37.bed sample1.bam > sample1_depth.txt
awk '{s+=$3}END{print s/NR}' sample1_depth.txt

There is already code in htslib for resolving overlapping paired-end reads, but it only seems to be used in samtools mpileup and not in samtools depth. I guess you could modify bam2depth.c and add the bam_mplp_init_overlaps() function call. I used it in my own program for strand-specific coverage calculation, see here

That is interesting indeed and I will definitely look into the bam2depth.c modifications in order to drop overlap from the function call. Was it difficult to implement this into your program? I see your program and the overlapping ends counting feature appears to be what I am looking for. I will need to locally install samtools as I work on a cluster. Thanks for the recommendation!

samtools depth is a simplified version of samtools mpileup, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use samtools mpileup instead:
$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...

I haven't shown the entirety of the output and you'd want to use the -aa flag to ensure that all bases are returned.
Note that samtools mpileup is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to -Q 1 (the default is -Q 13, which you'd want to change).
As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use bamPEFragmentSize from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.

Does this method take bases that are not covered by a read into account?

Yes, that's what the -aa flag I mentioned does.

I tested this on data and it seems to handle the issue nicely! Thanks Devon.

One way to deal with this would be to first merge paired-end reads based upon their overlapping regions, and then map them and calculate the coverage. This way you're only counting once per unique sequence.
Programs like SeqPrep, PEAR (Paired-End reAd mergeR), and fastq-join can do this fairly quickly. However, this isn't a perfect solution as the reads will not be joined unless they meet some minimum overlap criteria.

added an edit, I do not want to remap reads or make any modifications to existing reads.

picard CollectHsMetrics has the parameter PER_BASE_COVERAGE which take this into account.
You can run it like this:
 $ java -jar picard.jar CollectHsMetrics \
      I=input.bam \
      O=hs_metrics.txt \
      R=reference_sequence.fasta \
      BAIT_INTERVALS=bait.interval_list \
      TARGET_INTERVALS=target.interval_list \
      PER_BASE_COVERAGE=base_coverage.csv

BAIT_INTERVALS and TARGET_INTERVALS could be the same file. You will need to convert your bed file to it using picard BedToIntervalList.
Another useful tool is mosdepth.
$ mosdepth --by regions.bed SampleName input.bam

fin swimmer

Forgive me because I’be only looked into the tool, never used it. Do I have to generate some kind of sample.dict file for every sample I want to do this on or is it straightforward. It’s no big deal doing this on a reference but it would have too much overhead running it on every sample. Could you provide the steps involved in running this?

A sample.dict isn't needed. I edit my answer as @terdon asks for.

Bases from overlapping reads of the same fragment provide extra information, as their base qualities come from different sequencing cycles on the sequencer. This information can be used by downstream tools. For example, the merging tool of USEARCH increases base quality scores of bases that match between overlapping reads. Depending on the context under which coverage information will be used in the future, it may be a good idea to "double-count" overlapping bases.

The reads aren't from different clusters in the traditional sense, if that were the case they wouldn't be assigned as mates. The traditional method to exploit overlapping reads is to increase/decrease phred scores on one of the mates according to whether the base calls in overlapping regions are the same or different, respectively. That's performed internally in things like samtools mpileup.

I agree. The reads come from the same cluster.  However, you could match up mates from different clusters. The Illumina sequencer uses a fragment index to match up mates from the two mate sequencing runs. It also appears from your comment that samtools incorporates the extra information gained from having two reads cover the same base. Just so we're all on the same page, here is a video describing Illumina paired-end sequencing: https://youtu.be/fCd6B5HRaZ8


text:
Getting data from fastq by generator
I have a task in a training that I have to read and filter the 'good' reads of big fastq files.
I downsampled, got the code working, saving in a python dictionary. But turns out the original files are huge and I rewrite the code to give a generator. It did work for the down-sampled sample. But I was wondering if it's a good idea to get out all the data and filtering in a dictionary. Does anybody here a better idea? I am asking because I am learning python on my own and I am not sure if it makes sense.
I got some ideas from a code in Biostar:
import sys
import gzip

filename = sys.argv[1]

def parsing_fastq_files(filename):

    with gzip.open(filename, "rb") as infile:
        count_lines = 0
            for line in infile:
            line = line.decode()
            if count_lines % 4 == 0:
                ids = line[1:].strip()
                yield ids
            if count_lines == 1:
                reads = line.rstrip()
                yield reads
        count_lines += 1

total_reads = parsing_fastq_files(filename)
print(next(total_reads))
print(next(total_reads))

I now need to figure out to get the data filtered by using 'if value.endswith('expression'):' but if I use a dict for example, but that's my doubt because of the number of keys and values.

You want lines where the header endswith expression correct?

Have a function that returns a single-end read (or pair) and another that does any relevant filtering on that. Wrap both of them in another function that checks if you've hit the end of the file. That's both simple, will handle any size of file, and allows single and paired-end reads.

This is more of a comment than an answer right Devon?

@d_kennetz Not really, it answers what the proper scalable method is. One could include the actual code, but I have faith in the poster's ability given a prose description.

I do not have an example of fastqs with 'expression' as the last string in the header, so instead I will filter based off of a specific barcode used for demultiplexing. I have a gzipped fastq of 250 reads in size. Based on Devon's answer, I have first generated a fastq parser which stores only (header, seq) as a list of tuples, and then looks inside the list of tuples for a specific index. Some of the indexes vary by 1 nucleotide so those will be left out. I can prove this by printing the length of the 2 lists.
Here is the code:
import sys
import gzip

fastq = sys.argv[1]

def parser(filename):
    """
    Generate a list of tuples (header, read)
    """
    fastq_parsed = []
    try:

        with gzip.open(filename) as fq:
            header = next(fq)
            read = next(fq)
            fastq_parsed.append((header[1:-1], read[:-1])) # we don't want new-line chars or @ sign in reads
            while True:
                next(fq) # skip + line
                next(fq) # skip qscore
                header = next(fq) # store next read header
                read = next(fq) # store next read seq
                fastq_parsed.append((header[1:-1], read[:-1]))
    except:
        StopIteration # when we run out of reads, stop
    return fastq_parsed

def expression_finder(read_list):
    """
    Find headers that end with expression assuming the list format generated above.
    """
    filtered_reads = []
    for header, read in read_list:
        if (header.endswith(b'ATCACGAT')): #gzip files in bytes format, so we pass bytes object.
            filtered_reads.append((header.decode("utf-8"), read.decode("utf-8"))) # list items decoded to utf-8
    return filtered_reads

read_list = parser(fastq)
print(len(read_list))
subset_reads = expression_finder(read_list)
print(len(subset_reads))

Which returns:
(base) [dkennetz@nodecn203  test]$ python3.6 fastq_parser.py subset.fastq.gz
250
173

We can see that of the 250 reads, 173 have the correct index. I have hardcoded the filter in the second function, but you could also pass it as an argument to the function, like:
def expression_finder(read_list, filter):
    ...
        if header.endswith(filter):
            ...

subset_reads = expression_finder(read_list, b'expression')

which would also work. But include the b'' in your string syntax because the first function still passes a bytes object.

Hey d looks cool but if you .decode() the header and the read when you append? Thanks for share your ideas. paulo

I did decode when I append but you certainly don’t have to do so! You can leave them as bytes if you’d like by simply passing .append((header, seq)). This will leave each as bytes objects.


text:
How to separate a microbial strain of MTB resistant and susceptible drugwise?
Im trying to classify my list into three type: XDR - A strain must be resistant to:
MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN + AMYCACIN/CAPREOMYCIN/KANAMYCIN-> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN) + any of these or two or three or all (AMYCACIN/CAPREOMYCIN/KANAMYCIN)

PRE-XDR - A strain must be resistant to:
MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN -> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN)

MDR - A strain must be resistant to:
ISONIAZID + RIFAMPICIN
and if it is resistant to any one/two/three of these or all of these in
addition to ISONIAZID + RIFAMPICIN than again it is MDR:-
PYRAZINAMIDE, ETHIONAMIDE, STREPTOMYCIN, ETHAMBUTOL

SUSCEPTIBLE - A strain must be Susceptible to all drugs
I tried the code as follows for MDR:
tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P "\t,ISONIAZID,RIFAMPICIN$"

It gave me MDR starins correctly but than i had to first remove all Susceptible strains and applied it to Resistant one only.
But when i applied the code for Pre-XDR:
tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P
"\t,ISONIAZID,RIFAMPICIN|CIPROFLOXACIN|MOXIFLOXACIN|OFLOXACIN|LEVOFLOXACIN"

It is giving me the result as below:
SRR671821       ,CAPREOMYCIN,ETHAMBUTOL,ETHIONAMIDE,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671863       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN
SRR671822       ,ETHAMBUTOL,ISONIAZID,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671864       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN

Which is not correct for Pre-XDR.
So I'm looking for a the drugs if match than look for resistant and susceptible if these criteria match with any of the mention type (MDR,XDR...) it put those under that in a new file. Can anyone help?
my file druglist.sxls
MDR eg from the file attached is:
ERR038736   ISONIAZID   RESISTANT
ERR038736   RIFAMPICIN  RESISTANT
ERR038737   ISONIAZID   RESISTANT
ERR038737   RIFAMPICIN  RESISTANT
SRR1270146  AMIKACIN    SUSCEPTIBLE
SRR1270146  CAPREOMYCIN SUSCEPTIBLE
SRR1270146  CYCLOSERINE SUSCEPTIBLE
SRR1270146  ETHAMBUTOL  RESISTANT
SRR1270146  ETHIONAMIDE SUSCEPTIBLE
SRR1270146  ISONIAZID   RESISTANT
SRR1270146  KANAMYCIN   SUSCEPTIBLE
SRR1270146  PARA-AMINOSALICYLIC ACID    SUSCEPTIBLE
SRR1270146  PYRAZINAMIDE    RESISTANT
SRR1270146  RIFAMPICIN  RESISTANT
SRR1270146  STREPTOMYCIN    RESISTANT

Why did you delete your original question hat already had an answer?

Hi @terdon I didn't deleted that question and the answers were removed by itself and it went to solved category it self thats why i have to repost the question.

Well the answer was deleted by the user who posted it, yes, but the question was deleted by you.  Note that there is no such thing as "solved category", I am not sure what you mean by that. Anyway, don't worry about it, just try not to delete and repost questions in general since it confuses things. But no big deal.

I deleted my answer as it was incorrect (your question was correct and could receive answers from other people). I could try to answer it again with base R, but I will need more time

Please [edit] your question and show us an example of drug.txt

I have editted my question and have put a small data from my drug.csv file. I have also uploaded the whole file as well.

@llrs Thankyou. Any help would be appreciable.

Also asked 6 weeks ago on biostars: https://www.biostars.org/p/344050/

Since this involves combining multiple conditions, it is easier (and certainly more readable if you ever need to come back to it later) to write a script to do this. For example, in perl:
#!/usr/bin/perl

use strict;
use warnings;

## This will hold the information for each strain
my %strains;

## These are the functions that will test each strain
sub isMDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}){
    foreach my $drug (keys(%thisStrain)) {
      next if $drug eq "ISONIAZID";
      next if $drug eq "RIFAMPICIN";
      next if $drug eq "PYRAZINAMID";
      next if $drug eq "ETHAMBUTOL";
      next if $drug eq "STREPTOMYCIN";
      ## discard this strain if it is resistant to any other drug
      if ($thisStrain{$drug}) {
        return 0;
      }
    }

    return 1;
  }
  return 0;
}

sub isPreXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}) &&
  (
   $thisStrain{"CIPROFLOXACIN"}  ||
       $thisStrain{"MOXIFLOXACIN"} ||
   $thisStrain{"OFLOXACIN"}  ||
       $thisStrain{"LEVOFLOXACIN"}
      )
     ) {
    return 1;
  }
  return 0;
}

sub isXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  isPreXDR($_strain) &&
      (
       $thisStrain{"AMYCACIN"} ||
   $thisStrain{"CAPREOMYCIN"} ||
       $thisStrain{"KANAMYCIN"}
      )
     ){
    return 1;
  }
  return 0;
}

sub isSusceptible{
  my $_strain = shift;
  foreach my $drug (keys(%{$strains{$_strain}})) {
if ($strains{$_strain}{$drug}) {
      return 0
    }
  }
  return 1
}

## Iterate over each line in the input file.
## For each iteration, the special variable $_
# will hold the contents of the line.
while (<>) {
  ## Skip first line
  next if $. == 1;
  ## remove Windows-style line endings
  s/\r//;
  ## remove trailing newline
  chomp;
  ## split the line on commas
  my ($strain, $drug, $status) = split(/,/);
  if ($status eq "SUSCEPTIBLE"){
$status = 0
  }
  else {
$status = 1;
  }
  $strains{$strain}{$drug}=$status;
}

foreach my $strain (keys(%strains)) {
  print "$strain\t";
  if (isMDR($strain)) {
print "MDR\n";
  }
  elsif (isXDR($strain)) {
print "XDR\n"
  }
  elsif (isPreXDR($strain)) {
print "PRE-XDR\n"
  }
  elsif (isSusceptible($strain)) {
    print "SUSCEPTIBLE\n";
  }
  else {
    print "OTHER\n";
  }
}

If you save that script as foo.pl, you can run it on your file with:
perl foo.pl drug.csv

Selected output from running this on the drug.csv file you posted is:
$ foo.pl drug.csv | grep -E 'ERR038288|ERR038736|ERR038751|SRS553990|ERR133879'
ERR038751   OTHER
ERR038288   SUSCEPTIBLE
SRS553990   XDR
ERR038736   MDR
ERR133879   PRE-XDR

Various details about this question were added in chat: https://chat.stackexchange.com/transcript/message/47821646#47821646

It is again not correct like for MDR -ERR047000 - It is a Susceptible strain.
ERR047000 ETHAMBUTOL SUSCEPTIBLE
ERR047000 ISONIAZID SUSCEPTIBLE
ERR047000 PYRAZINAMIDE SUSCEPTIBLE
ERR047000 RIFAMPICIN SUSCEPTIBLE
Same for MDR -> ERS389515:
ERS389515 ETHAMBUTOL SUSCEPTIBLE
ERS389515 ISONIAZID SUSCEPTIBLE
ERS389515 PYRAZINAMIDE SUSCEPTIBLE
ERS389515 RIFAMPICIN SUSCEPTIBLE
ERS389515 STREPTOMYCIN SUSCEPTIBLE
It is susceptible not MDR

@SafinaA.R what is wrong? It is very hard to test since you haven't given examples of what we should see. The example in your question isn't listed in the drugs.csv file. Can you ping me (@terdon) in our chat room so we can discuss this?

HI. I think it is now running fine. So now i have the proper classification verified by my team. So any thing that is resistant to ISONIASID and RIfampicin is mdr but it can also be mdr if these two are there plus any of these or all of these: Pyrazinamid, Ethambutal, Streptomycin, ETHIONAMIDE. like these are the few example of mdr from the list:
ERR038736 ERR038737 ERR038738 ERR038739 ERR038740 ERR038742 ERR038743 ERR038744 ERR038745 ERR038746 ERR038747 ERR038748 ERR038749 ERR038750 ERR038752 ERR046855

Few Examples:
ERR038288 Susc ERR038736 MDR ERR038751 Others SRS553990 XDR ERR133879 Pre-XDR

@SafinaA.R all of those are caught correctly by the latest version.

Great. Correctly classified. THUMBS UP.


text:
summarising read group information from a .bam file
I have merged together 2 different .bam files in order to simulate sample contamination. So the reads can come from one of two samples, as shown by the read group info:
@RG ID:0    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:1:none
@RG ID:1    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:2:none
@RG ID:2    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:3:none
@RG ID:3    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:4:none
@RG ID:0-11EFC00B   PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:1:none
@RG ID:1-B8A1099    PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:2:none
@RG ID:2-330086F    PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:3:none
@RG ID:3-7681F092   PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:4:none

I'd like to check that the correct proportion of read groups originate from each sample.
Currently I am using:
samtools view example.bam | rev | cut -f 1 | rev > output.txt

, but this is not very elegant and only works because the RG field is last in the .bam.
Is there a quick way to tabulate the number of reads groups with different IDs? E.g. produce an output like:
ID:0 1000
ID:1 2000
ID:2 3000
...

A solution in samtools would be ideal, along the lines of the output produced in samtools stats.

Does this answer your question: https://bioinformatics.stackexchange.com/q/10559/6251 ?

Also, try sambamba.

@TimurShtatland, yes, using samtools view -R read_groups.txt example.bam | wc -l will work. thanks.

Would you like to self-answer the question then? So it's easier for others to find the solution.


text:
Annotating a .vcf with centimorgan information
Some programs (e.g. shapeit4) automatically annotate an INFO tag into a .vcf file which gives the cumulative genetic distance in cM between each SNP:
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=10/11/2020 - 21:29:39
##source=G2H
##contig=<ID=9>
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=1,Type=Integer,Description="Allele count">
##INFO=<ID=CM,Number=A,Type=Float,Description="Interpolated cM position">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Phased genotypes">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##bcftools_viewVersion=1.9-207-g2299ab6+htslib-1.9-271-g6738132
##bcftools_viewCommand=view -s -1_-1 ukb_cal_chr9_v2.HumanOrigins.common_snps.HumanOriginsInds.fiftyPercentEthnicUKBInds.PHASED.bgz.vcf.gz; Date=Mon May 24 14:15:51 2021
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  -1_-1
9   257034  rs573860    A   G   .   .   AC=0;AF=0.388699;CM=0;AN=2  GT  0|0
9   270321  rs17720310  A   G   .   .   AC=0;AF=0.0820595;CM=0.007994;AN=2  GT  0|0

However, this program is really used for phasing not for .vcf annotation.
Ideally I'd like a simple command line program which takes input the .vcf and a genetic map and interpolates the genetic distance and adds a cumulative map. This could be done manually with a custom script, but quite frankly I can't be bothered to do that.
Does anyone know of a tool that exists which will do this?
I posted this to biostars a while ago and received one reply which linked to this post, which includes links to tables containing cM locations for different rs numbers.
EDIT: I am interested in doing this for human data (although the same would be the case for any diploid organism).
Sample vcf:
https://pastebin.com/raw/YSbuU0qM

I think a useful tool to address your question is QCTOOL v2. Particularly, the option you are looking for is -annotate-genetic-map. This tool has been developed originally for GWAS datasets. It has multiple inputs data format.
If you could provide input data as an example, I will try working on it to provide the necessary steps to be done.
Hope this answers your question.

Sorry took a while, have added the sample vcf

@user438383 could you please specify whether your VCF is of hg38 or hg19? Also, do provide a proper VCF format (it is missing some columns like FORMAT followed by another column like in your original post example ). I tried using the file you provided, it is giving an error of segmentation fault. When I tried to add a vcf header with the proper format aforementioned, it works.


text:
Why is the p-value significance threshold for HLA association tests $5*10^{-8}$?
Typically, a p-value of $p<=5*10^{-8}$ is used for genome-wide association testing, as there are roughly $1/p$ independent LD blocks in the human genome, so it correct for multiple testing across the genome.
HLA assocation tests, such as in this study, only look at a fraction of the genome, something like chr6:28–34 Mb, which is approximately 0.2% of the genome used to conduct a true 'genome-wide' association study. However, they still use the same genome-wide p-value threshold of $p<=5*10^{-8}$.
Is this not too conservative, since there is a much reduced multiple testing burden? Should the threshold not be something like $p/0.002$?

Its a very good question and it is also a good study. What you are saying must be correct.
I thought it was possible they used the correct critical probability in the original submission and a reviewer requested it was changed. However I checked their medRxiv upload here and the authors clearly use this threshold prior their formal submission.

We performed stepwise conditional analysis 479 to identify additional
independent signals by adjusting for the most significant amino acid
480 position in each step until none met the significance threshold (
P = 5 10^-8)

Thus I don't know.
If this paper is central to your area of interest I would use the corrected critical probability threshold and re-assess and extend their conclusions because 5 decimal places is a huge a difference.

As a direct answer to your question, I would argue that the authors, in choosing this threshold, have demonstrated a lack of understanding of their process of analysis. But that answer is not a full answer, and it's a problematic answer because [almost] everyone who does population genetics does similar things. This p-value is not appropriate, not merely because it is at a genome-wide significance level, but because no threshold is appropriate.
To explain my perspective more, thresholding to identify "statistically-significant" results, in general, is not a good approach. Regardless of what threshold is used, p-values only indicate the confidence in a statistic's value (as it relates to fitting a particular statistical model), not importance or the size of the value.
The authors of this paper, as with many other authors who have written similar papers, make an obvious statistical error in the legend of figure 3, where they state: "One-field classical allele HLA-B*57 (P = 9.84 × 10−138) (a) and amino acid position 97 in HLA-B (Pomnibus = 1.86 × 10−184) (b) showed the strongest association signal." This is an error, because they have only reported here the p-values, and p-values have no relevance to the strength or magnitude of a statistic.
Given this, it doesn't make sense to rank by p-value, and it especially doesn't make sense to display p-values on a graph and highlight p-value peaks as indicating association peaks. Those regions are confidence peaks.
In an attempt to provide some constructive light around this, I will answer the obvious follow-up questions are "How should p-values be used?" and "If p-values shouldn't be displayed or ranked, then what should?"
My preference is to use p-values as a loose filter, to exclude obviously incorrect results. In the case of a GWAS, I might, for example, choose to exclude any statistic that has a p-value above 0.1, leaving the remainder as potentially correct. It could possibly be argued that this is still thresholding. It is, in a sense, but importantly I am not saying that all values that pass the threshold are good, just that the ones that are excluded are all bad; I allow for some rubbish to end up in the good pile. I liken it to using a computer to find the things that look like diamonds, so that I don't waste my time inspecting dirt.
After that filtering is done, the p-values should be ignored, and the descriptive statistic used to rank, determine significance, association, or whatever the goal of the analysis is. For GWAS, a commonly-used value is beta, which can be positive or negative depending on the direction of association, and relates to the variant difference between groups.
If a list needs to be cut down, then this descriptive statistic should be used for thresholding, depending on the analysis. For example, it might be appropriate to display the 100 most extreme beta values, or threshold to exclude any beta values that have a magnitude less than 0.01, or display a manhattan-like plot with beta on the $y$ axis, and chromosome position on the $x$ axis.
There shouldn't be set rules, because every situation is different - there's a human element to this. What I'm saying is that it should be the job of the researcher, not the computer, to identify which data are significant. Make the data set small enough that it can be processed by a person, then let them do the rest. I think Hannah Fry explains this quite well in this talk:

Now this is the stuff the algorithms are amazing at: looking for tiny, tiny clues in seemingly completely disconnected data sets.... Just incredibly, incredibly sensitive these algorithms, right? Just absolutely amazing what they can do. And yet remember that sensitivity does not make a perfect algorithm. That's only one half of the equation.... [Humans] are terrible at sensitivity, right? We are really, really bad. We miss things all the time. But, specificity, being specific, that's like our superpower.... And rather than choosing between human and machine, which is kind of the rhetoric that we get so much, why don't we exploit each other's strengths and just create much more of a partnership?.... Because the algorithm never gets tired, so let it trawl through all of that data and just highlight a few key areas of concern. And then the human never misdiagnoses. So they can come in, and just sweep up and have the final say.


text:
Low custom tdtomato gene content
I have a set of scRNA-seq samples expressing TdTomato, which has high content in microscope. I followed the 10x cellranger pipelines to finsh the work, my procedures are as follows:

added TdTomato on an additional chromosome:

added an entry in the gtf
cellranger mkref
cellranger count

I got the result, which includes:

But strangely, the tdtomato expression content in my three sample are all very low! And cells with tdtomato are also very low and disperse!! Take sample 1 as example, you can see the picture:

It's not normal, but I can't find the reason. Following the 10x genomics help documentation noticed that by executing the CellRanger filtering steps, maybe my tdtomato reads are filtered too much, but because I'm not very skillful at programming and I'm first time to do RNA-seq (and scRNA-seq), I just don't know how to analysis the bam file as suggested by the official website and find the real reason of the problem:

If you can show me the code need to do, I will be very grateful!

Just because you can see the protein with the microscope doesn't mean that the RNA is expressed at a sufficient quantity to detect it in a single cell experiment.  Single cell RNA-seq generally skims the top most expressed genes.

but its rna level is obviously lower than average,and cells expressing it is few,this is totally wrong,so the align must have some wrong problems to solve.

Have you looked in the bam file to see if there are a lot fo reads aligning to tdTomato? I am also not entirely sure what you are saying.  I know of other people who have done a similar experiment and were unable to detect the gene because it's expression level was too low.

Possible duplicate of No counts for added gene in cellranger (scRNA-seq)

Any update on whether including the 3' UTR fixed the issue? Thank you!

A few suggestions:

Manually inspecting the BAM file with a genome browser like IGV, as @GWW suggested and checking the qualities of the alignments in your BAM/SAM file (SAM is the human readable version of BAM, there are tools to convert one to another);
Using another single cell analysis tool like Seurat so that you can manually set thresholds for QC and not lose to many (tdtomato) expressing cells at the QC step (I don't know about your experiment but if transfected, cells could be at a worse condition than non-transfected ones and hence might not pass QC);
Double checking your tdtomato reference sequence, I don't know about tdtomato but for example GFP has flavors like EGFP, similar but different sequences, might explain low alignment score if you observe in step 1;
Taking the limitations of scRNA-seq (especially droplet-based) technologies into account. The figure below is from the CITE-seq paper (from supp info), Simultaneous epitope and transcriptome measurement in single cells (Stoeckius et al., 2017), showing that the protein - RNA measurements of a given gene might show poor correlation.

The solution to your problem probably is adding the full mRNA sequence of your transgene to the reference (as also suggested by acrux).
I had a very similar problem recently when tdTomato expression was only detectable in a couple of cells. The transgene was introduced into the cells using a lentiviral vector and therefore the mRNA had the tdTomato sequence followed by a WPRE element. In the picture below you can see that a large number of reads is mapping to the WPRE element in the 3' UTR (blue region) and and almost none to the tdTomato CDS (red region).

In general, the untranslated regions are considered exons because they are part of the final mRNA. So if you have the fasta sequence beginning at the transcription start and ending at the polyadenylation signal, this is the full mRNA. Then you just need to update the GTF with the longer full mRNA length and you are all set.

PPK and and acrux, has any of this adding of the 3'UTR sequence to the coding for sequencing transgenes been published?

As I mentioned, the 3' UTR is part of all gene annotations. Without these sequences your reference is incomplete leading to the failure to detect the transgene.

We had the same issue. Turned out, our tdTomato transgene was followed up by WPRE and BGH polyA signals. These sequences are quite large, and their presence precludes sequencing into the tdTomato coding sequence when using the 10x genomics platform. Incorporating these into our build, allowed us to detect tdTomato+ cells.

I am unable to respond to other's comments as I never usually comment, so this is a new account. However, I spent the last couple days really diving into the BAM alignments when mapping scRNAseq reads to the human genome, but altering the GTF and FASTA using different annotations of my CAR construct (some very liberal allowing upstream and to the poly A tail sequence, and some less liberal)
MY MAIN CONCLUSION!:
DO NOT ADD THE BGH POLY A SIGNAL TO THE GENE ANNOTATION! it is highly homologous to the human beta globin subunit, and I am finding almost none of the reads actually mapping both to the adjacent CAR construct and the BGH POLY A subunit.
I am still figuring out how to handle the CD3zeta, CD28, and the OX40 in the CAR construct being multi-mapped (I will test just using a small portion of the CD28 transmembrane domain and then using just the chimeric antigen portion for mapping)
BUT, BIGEST TAKEAWAY, if you are going to attempt to add the poly A tail or any other genomic annotation not strictly relating to your gene of interest look at the bam and ensure that the reads mapped actually bridges into your gene. (as I also had residual genomic DNA from human beta globin upstream of my CAR construct that showed splicing over the CAR construct and aligning downstream of the CAR to the POLY A tail sequence)
-for viewing the BAM i am using IGVtools, and for the creation of the GTF and the FASTAs that I concatenated to the human reference genome I used tears [not a software package] ;)

you might want to try adding the 3'UTR sequence of you tdtomato construct used in your cells. that should do the trick.

Welcome to the site acrux. Thanks for answering and helping out. Could you expand a bit how does adding this helps to the original poster?


text:
length of 'dimnames' [1686] must match that of 'dims' [3]
Please if anyone has experience with the use of the BSEQ-SC package for the deconvolution of bulk RNA sequencing data with single cell RNA sequencing data I will be very grateful for your suggestion.
I wish to deconvolve my bulk RNA seq data with single cell RNA seq data from the same tissue. I run in to the following error when running the code line below and I do understand the error but not how I can overcome it. Can any one please provide any useful hints?
csfit <- bseqsc_csdiff(ilc_bulk[genes, ] ~ treatment | fibroblast + macrophage,
                       verbose = 2, nperms = 100, .rng = 12345)

Then comes the error:
   Groups: a_nacl1=4L | a_nacl2=4L | a_nacl3=4L | a_nacl4=4L | a_nacl5=4L | NA0L
Cell type(s): 'fibroblast', 'macrophage' (2 total)
Fitting mode: auto
Data (filtered): 1684 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ...  OK
 Computing FDR using 100 permutations ... 101/100
  Alternative 'two.sided' ...   OK
  Alternative 'greater' ...   OK
  Alternative 'less' ...   OK
 OK
Timing:
   user  system elapsed
  2.136   0.153   2.267
Error in array(NA_real_, dim = c(nrow(b), ncol(cc), length(lev)), dimnames = c(rownames(b),  :
  length of 'dimnames' [1687] must match that of 'dims' [3]

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] edgeR_3.26.6                limma_3.40.6                forcats_0.4.0               stringr_1.4.0
 [5] purrr_0.3.2                 readr_1.3.1                 tidyr_0.8.3                 tibble_2.1.3
 [9] tidyverse_1.2.1             SoupX_0.3.1                 ggplot2_3.2.0               dplyr_0.8.3
[13] DropletUtils_1.4.3          SingleCellExperiment_1.6.0  SummarizedExperiment_1.14.0 DelayedArray_0.10.0
[17] BiocParallel_1.18.0         matrixStats_0.54.0          GenomicRanges_1.36.0        GenomeInfoDb_1.20.0
[21] RColorBrewer_1.1-2          xbioc_0.1.17                AnnotationDbi_1.46.0        IRanges_2.18.1
[25] S4Vectors_0.22.0            BisqueRNA_1.0               Seurat_3.0.2                preprocessCore_1.46.0
[29] e1071_1.7-2                 bseqsc_1.0                  csSAM_1.4                   Rcpp_1.0.2
[33] openxlsx_4.1.0.1            Biobase_2.44.0              BiocGenerics_0.30.0

loaded via a namespace (and not attached):
  [1] reticulate_1.13        R.utils_2.9.0          tidyselect_0.2.5       RSQLite_2.1.2          htmlwidgets_1.3
  [6] grid_3.6.1             Rtsne_0.15             devtools_2.1.0         munsell_0.5.0          codetools_0.2-16
 [11] ica_1.0-2              future_1.14.0          withr_2.1.2            colorspace_1.4-1       rstudioapi_0.10
 [16] ROCR_1.0-7             gbRd_0.4-11            listenv_0.7.0          NMF_0.22               Rdpack_0.11-0
 [21] labeling_0.3           GenomeInfoDbData_1.2.1 bit64_0.9-7            rhdf5_2.28.0           rprojroot_1.3-2
 [26] vctrs_0.2.0            generics_0.0.2         R6_2.4.0               doParallel_1.0.14      rsvd_1.0.2
 [31] locfit_1.5-9.1         bitops_1.0-6           assertthat_0.2.1       SDMTools_1.1-221.1     scales_1.0.0
 [36] gtable_0.3.0           npsurv_0.4-0           globals_0.12.4         processx_3.4.1         rlang_0.4.0
 [41] zeallot_0.1.0          splines_3.6.1          lazyeval_0.2.2         broom_0.5.2            modelr_0.1.4
 [46] BiocManager_1.30.4     yaml_2.2.0             reshape2_1.4.3         backports_1.1.4        tools_3.6.1
 [51] usethis_1.5.1          gridBase_0.4-7         gplots_3.0.1.1         sessioninfo_1.1.1      ggridges_0.5.1
 [56] plyr_1.8.4             zlibbioc_1.30.0        RCurl_1.95-4.12        ps_1.3.0               prettyunits_1.0.2
 [61] pbapply_1.4-1          viridis_0.5.1          cowplot_1.0.0          zoo_1.8-6              haven_2.1.1
 [66] ggrepel_0.8.1          cluster_2.1.0          fs_1.3.1               magrittr_1.5           data.table_1.12.2
 [71] lmtest_0.9-37          RANN_2.6.1             fitdistrplus_1.0-14    pkgload_1.0.2          hms_0.5.0
 [76] lsei_1.2-0             xtable_1.8-4           readxl_1.3.1           gridExtra_2.3          testthat_2.2.1
 [81] compiler_3.6.1         KernSmooth_2.23-15     crayon_1.3.4           R.oo_1.22.0            htmltools_0.3.6
 [86] Formula_1.2-3          lubridate_1.7.4        DBI_1.0.0              MASS_7.3-51.4          Matrix_1.2-17
 [91] cli_1.1.0              R.methodsS3_1.7.1      gdata_2.18.0           metap_1.1              igraph_1.2.4.1
 [96] pkgconfig_2.0.2        registry_0.5-1         plotly_4.9.0           xml2_1.2.1             foreach_1.4.7
[101] dqrng_0.2.1            rngtools_1.4           pkgmaker_0.28          XVector_0.24.0         rvest_0.3.4
[106] bibtex_0.4.2           callr_3.3.1            digest_0.6.20          sctransform_0.2.0      tsne_0.1-3
[111] cellranger_1.1.0       dendextend_1.12.0      curl_4.0               gtools_3.8.1           nlme_3.1-141
[116] jsonlite_1.6           Rhdf5lib_1.6.0         desc_1.2.0             viridisLite_0.3.0      pillar_1.4.2
[121] lattice_0.20-38        httr_1.4.0             pkgbuild_1.0.3         survival_2.44-1.1      glue_1.3.1
[126] remotes_2.1.0          zip_2.0.3              png_0.1-7              iterators_1.0.12       bit_1.1-14
[131] class_7.3-15           stringi_1.4.3          HDF5Array_1.12.1       blob_1.2.0             caTools_1.17.1.2
[136] memoise_1.1.0          irlba_2.3.3            future.apply_1.3.0     ape_5.3

Answer from @stupidwolf, @benn, @llrs, converted from comments:
The error complains about the 1686 rownames, while you have 1684 features. I think something is wrong with the naming of your rows.
This may be a bug that you should report to the package maintainer. I cannot really tell without more information about your eset. Looking at this line "Cell type(s): 'fibroblast' (1 total)", is it possible you only have one group under fibroblast? It will not make sense to run a random effect for only 1 group.
[please edit to improve this answer, if possible]


text:
Gene expression changes in one condition
I want to look at an RNA-seq data-set on 1 condition to look for expressed genes. As there is only one condition, I can’t perform differential analysis.
The only thing I can think to do is set a threshold (i.e. all replicates must have > 10 gene counts) then say these are expressed. But I don’t think this is a very good way of doing it.
Do you know of any better ways that we can do this?

If you've got enough replicates (ideally at least 6, but can be done with 3), you should be able to do a Variance-stabilised Transformation using DESeq2:
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
head(assay(vsd), 3)

More details on this can be found in the DESeq2 Vignette:
https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#extracting-transformed-values

Answer from @james-hawley, converted from comment:
In any high throughput sequencing experiment, there will be lowly expressed RNA molecules that are not captured by the library prep + sequencing. So defining what "is expressed" is somewhat arbitrary from RNA-seq data, regardless of how many replicates you have and what method you use. That said, you can control some of the noise in your RNA-seq data using a variance-stabilized transformation, like @gringer suggests, and look at the distribution of genes/transcripts. Some threshold like $rlog \ge x$ in at least $n$ samples is usually a good starting point (for some arbitrary $x$ and $n$).

I don't think there is a very good way to do what you want to do.  I'd think some kind of RPM would be better than strictly raw counts; you have to account for overall library size.  Or TPM might be better.


text:
Unifying these files together
I have two .csv files
One for raw read counts of each gene and the other 3114 single cells for four patients which should be the row names of the first file
As single cells per patients seems duplicated I am getting error and don't know how to read them together
Any way this is the files
https://www.dropbox.com/s/slkhspz7lgdgph6/ascites.csv?dl=0
https://www.dropbox.com/s/iq8lt7qrh7hdtq8/ascitesRownames.csv?dl=0
I tried
a=read.csv("ascites.csv")
Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  restarting interrupted promise evaluation
b=read.csv("ascitesRownames.csv")
View(b)
View(a)
rownames(a)=b

Error in .rowNamesDF<-(x, value = value) :    duplicate 'row.names'
are not allowed In addition: Warning message: non-unique value when
setting 'row.names': ‘logical(0)’

Any help in having a matrix in row names are 3114 single cells and genes in columns?

Can you please dput() the head of files so we don't have to download them from dropbox? Thanks.

There are a couple of problems:

As Ram RS said you get a data frame when you load that second CSV.  Since the information is stored in a row rather than a column you need a little more work to extract it too.
R will silently mangle your column names from the first data frame and will also interpret the second file as having column names (and then mangle those too) unless you give extra arguments to read.csv to prevent it.
The fact that you only have four identifiers (not 3114) in the second file means you can't use those as row names.  Your rows in a are cells, not patients, right?  You could make up something unique based on both patient and cell info but I'd suggest avoiding that.

First you could use something like this (going straight from Dropbox here) to read the files without losing information:
a <- read.csv("http://dl.dropboxusercontent.com/s/slkhspz7lgdgph6/ascites.csv", check.names = FALSE)
b <- read.csv("http://dl.dropboxusercontent.com/s/iq8lt7qrh7hdtq8/ascitesRownames.csv", header = FALSE)

Then, since you want to associate more information than can easily fit in a "wide-format" (matrix-like) data frame, and row and column names are so easily mangled by R, I'd suggest using a "long-format" data frame instead.  The package reshape2 is very helpful for this.  Here you can put the patient IDs as an additional column in your first data frame, explicitly label the cells, and then make a new data frame where each row is a cell and gene combination for a particular patient:
a <- cbind(patient = as.character(b[1, ]), cell_id = 1:nrow(a), a)
gene_counts <- reshape2::melt(
  a, id.vars = 1:2, variable.name = "gene", value.name = "count")

So you get something like this:
> head(gene_counts)
  patient cell_id gene count
1   7873M       1 A1BG     1
2   7873M       2 A1BG     0
3   7873M       3 A1BG     2
4   7873M       4 A1BG     0
5   7873M       5 A1BG     1
6   7873M       6 A1BG     0

(Note the gene column is a factor there.)
I'd be very careful with row and column names since R is so prone to altering them without warning, and in general with assuming things will match up in the right order like the rows of the one table and the columns of the other.  Hadley Wickham has a good summary of the problems with row names in particular:
https://adv-r.hadley.nz/vectors-chap.html#rownames

+1 agree with the comment about row names. I think data.table has the right idea in doing away with them completely and using keys instead.

read.csv produces a data.frame. You'll need rownames(a) <- b$col_name, not rownames(a) <- b.


text:
Parsing .vcf file for this information
I have a .vcf file
https://www.dropbox.com/s/8v73nppwg3a1tnd/LP2000109-DNA_A01_vs_LP2000103-DNA_A01.SVannotated.txt?dl=0
with this header
##startTime=Fri Mar 29 16:46:32 2019
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  NORMAL  TUMOR
1   54586   .   T   C   .   PASS    DP=39;MQ=50.55;MQ0=0;NT=ref;QSS=48;QSS_NT=48;ReadPosRankSum=1.92;SGT=TT->CT;SNVSB=0.00;SOMATIC;SomaticEVS=10.83;TQSS=1;TQSS_NT=1    AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,0:0,0:20:0:0,0:0:0:20,20  0,0:6,6:18:0:0,0:0:0:12,13
1   103241  .   C   T   .   PASS    DP=120;MQ=24.94;MQ0=35;NT=ref;QSS=47;QSS_NT=47;ReadPosRankSum=2.09;SGT=CC->CT;SNVSB=0.00;SOMATIC;SomaticEVS=9.44;TQSS=2;TQSS_NT=2   AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,1:32,47:33:1:0,0:0:0:0,5  0,

The "DP" field in the vcf shows the depth of the individual samples; So in this file, the first locus has the following format fields:
AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,0:0,0:20:0:0,0:0:0:20,20  0,0:6,6:18:0:0,0:0:0:12,13

So according to this (DP field of normal and tumor samples), normal sample has a depth of 20 and tumor sample has a depth of 18.
So how I could extract the read depth for all loci as described for the first position? The desired output would be like this [note that the VCF is taken from my own data but the table is my desired format that I don't know how to get that from my own data. chr have been added manually because my reference genome is hg19]:
  Sample   Type   CHROM       POS    REF  ALT       Tumor_Depth Normal_Depth
  CHC2432T SNV    chr1  102961055    G    A                  64           62
  CHC2432T SNV    chr1  105492588    A    T                  66           73
  CHC2432T SNV    chr1  108628724    C    T                  45           54
  CHC2432T SNV    chr1  109692113    G    T                  53           29
  CHC2432T SNV    chr1  109692114    G    T                  53           31
  CHC2432T SNV    chr1  120676701    T    C                  48           87

Your vcf file is truncated (see tail -n1 _double_filtered_013_pre_indels.vcf). I assume that's just the copy you uploaded to dropbox, but check your actual file to be sure.

You can do the extraction part with the GATK tool VariantsToTable, as described here:
https://gatk.broadinstitute.org/hc/en-us/articles/360041414592-VariantsToTable
The usage example from that doc:
gatk VariantsToTable \
    -V input.vcf \
    -F CHROM -F POS -F TYPE -GF AD \
    -O output.table

would produce a file that looks like:
 CHROM  POS        TYPE   HSCX1010N.AD  HSCX1010T.AD
 1      31782997   SNP    77,0          53,4
 1      40125052   SNP    97,0          92,7
 1      65068538   SNP    49,0          35,4
 1      111146235  SNP    69,1          77,4

So you might still need to reorder columns etc but that should allow you to at least get the values out in a tabular format that will be easier to work with.

To extract the DP fields from a VCF file, you could use a tool like bcftools query:

Extracts fields from VCF or BCF files and outputs them in user-defined format.

You could start from something like this:
bcftools query -Hf 'CHC2432T\t%TYPE\t%CHROM\t%POS\t%REF\t%ALT[\t%DP]\n' file.vcf

Sorry, how I can add RU=TAC part to tumour sample?


text:
Filtering a vcf with a text file of SNP rsIDs
I have a vcf file containing an ID field with a snps and a list of human SNPs from dbSNP.I want to remove whatever in dbSNP from vcf file. The text file is like this:
rs978760828
rs1008829651
rs1052373574
rs1326880612
rs768019142
rs892501864
rs1010989343
rs1489251879
rs1022805358
rs1322538365

How to remove the common snps from dbSNP and my vcf?

I'd advise against using grep, sed, awk for processing vcfs unless you really know what you're doing - there's almost always an easier way using bcftools.
bcftools view -i'ID=@snplist.txt' in.vcf


text:
Error in DESeq2 analysis with dimnames
I'm getting the following error in a DESeq2 analysis and I can't seem to figure out the issue:
Error in `assays<-`(`*tmp*`, withDimnames = withDimnames, ..., value = `*vtmp*`) :  please use 'assay(x, withDimnames=FALSE)) <- value' or 'assays(x, withDimnames=FALSE)) <- value' when the dimnames on the supplied assay(s) are not identical to the dimnames on RangedSummarizedExperiment object 'x'
Here is my code:
library(DESeq2)
library(edgeR)

url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
           "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
           "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
           "GSM1545545_JMS9-P8c.txt")

x <- readDGE(files, columns=c(1,3))

samplenames <- substring(colnames(x), 12, nchar(colnames(x)))

colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
                     "Basal", "ML", "LP"))
x$samples$group <- group

#DESeq2 Analysis
phenodata <- new("AnnotatedDataFrame", data = data.frame(celltype = x$samples$group))
counts <- x$counts
ids <- c(1,3,4,6,7,9)

tmpcounts = counts[, ids]
tmpph = as(phenodata[ids, ], "data.frame")
rownames(tmpph) <- colnames(tmpcounts)

dds = DESeqDataSetFromMatrix(countData = tmpcounts, colData = tmpph, design = ~ tmpph$celltype)

Here are my package versions
> packageVersion("DESeq2")
[1] ‘1.30.1’

> packageVersion("SummarizedExperiment")
[1] ‘1.20.0’


text:
Why index a sequence file with minimap2?
I have been using minimap2 for a few month now and I just found a script that follows these steps:
minimap2 [-x preset] -d target.mmi target.fa
minimap2 -a [-x preset] target.mmi query.fa > output.sam

I have never seen this .mmi index before. I always used minimap2 directly like this:
minimap2 -a [-x preset] target.fa query.fastq > output.sam

So my question is why should I index my reference genome before mapping? Is it only for speed up? Will it up speed my code if I do it? Or it is already done without specifying it?

If you don't provide an index, minimap2 will just generate the index on the fly (so, yes, it's done without specifying it). As such the mapping step will take longer. However, index generation is supposed to be very fast in minimap2 (a few minutes for the human genome), so it probably won't matter much.
Doing it separately will therefore also not speed up your code unless you end up using the pre-generated index multiple times.


text:
VCF file inspector app with GUI
Does anyone know a good, free-to-use VCF file inspector app (Mac OS, with GUI) where I can browse the content of a VCF file similarly to bamseek*?
I would like to browse SNP annotation split by columns, possibly even filter it.
'* Bamseek does not work for me, fails at open file. It anyways seems to be an abandonware.
Example table view from an expensive commercial software:

I wrote vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html
bcftools view in.bcf |  java -jar dist/vcf2table.jar | more

and http://lindenb.github.io/jvarkit/SwingVcfView.html

java -jar dist/swingvcfview.jar indexed.vcf.gz

OK it's a command line program, not a GUI, but bcftools query and bcftools view can easily be used to to view elements of a vcf/bcf file and are (probably) a lot more flexible and scalable than a GUI. Examples from here:
$ bcftools query -f '%CHROM %POS %REF %ALT\n' file.bcf | head -3
1 13380 C G
1 16071 G A
1 16141 C T

$ bcftools query -f '%CHROM %POS  GTs:[ %GT]\t PLs:[ %PL]\n' file.bcf | head -3
1 10234  GTs: 1/1 1/1    PLs: 28,3,0 29,3,0
1 10291  GTs: ./. 1/1    PLs: 0,0,0 28,3,0
1 14907  GTs: 0/1 0/1    PLs: 8,0,17 26,0,48


text:
Is it possible to cluster a sc-RNAseq expression matrix with weighting or embedding?
I have an RNAseq mouse retina gene expression matrix(24659,49300) to cluster effectively. Here is a sample(7,4):
gene            cell1   cell2   cell3
KITL            0       0       1
TMTC3           3       0       0
CEP290          1       3       0
4930430F08RIK   2       1       2
1700017N19RIK   0       0       0
MGAT4C          0       0       4

Is there any chance to apply a weighting schema or embedding to increase ARI, NMI or DBI? Dataset is Macosko, 2015: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63472
I want to try a new approach to cluster rnase data. Question is effectiveness to use weighting schema or embedding. Is it worth to try? May be someone tried it before.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.


text:
Combine VCF files
I have a problem with using rbind to combine VCF files using the library VariantAnnotation from Bioconductor.
I am reading two VCF files, when I try to combine them in a certain order with rbind I'm getting an error. When I combine them in a different order, it works fine.
vcf.full <- rbind(vcf1, vcf2) ## Doesnt work

Error in as(from, to_class, strict = FALSE) :    no method or default
for coercing “CompressedCharacterList” to “DNAStringSetList”

vcf.full <- rbind(vcf2, vcf1) ## Work

header(vcf1)
class: VCFHeader
samples(503): NA06984 NA06989 ... NA20828 NA20832
meta(2): META contig
fixed(2): FILTER ALT
info(27): CIEND CIPOS ... EX_TARGET MULTI_ALLELIC
geno(1): GT

header(vcf2)
class: VCFHeader
samples(503): NA06984 NA06989 ... NA20828 NA20832
meta(2): META contig
fixed(2): FILTER ALT
info(27): CIEND CIPOS ... EX_TARGET MULTI_ALLELIC
geno(1): GT

dim(vcf1)
[1]  30 503

dim(vcf2)
[1] 149 503

The two VCF files seem identical. Any idea of why it happens? Or any suggestions on a different method to combine multiple VCF files using R?

Answer from @Kozolovska, converted from comment:
Found the problem, some of ALT were structural variant, so the list of vcf2@info$ALT was a charterer list and not a 'DNAStringSetList' after removing that variant it worked. Thanks.


text:
How can I get latest .vcf files with annotation data?
I'd like to perform the annotation (1000genomes, COSMIC etc.) of my variants using SnpSift and SnpEff, however so far all I get are vcf files for separate chromosomes:
http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
I'd like to annotate my variants and then filter out common SNP's. What in your opinion is best annotation source and where can I get the file?

Do you mean Somatic Sequence Variants?,
if so, Table 1 inside the following full-text link:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5707196/
And additional information from my search-->
Here's
https://docs.google.com/document/d/1hVtgMYAtvDoUcigKx8OLDmongDRzDQAd/edit?usp=sharing&ouid=102408774767674945852&rtpof=true&sd=true

https://drive.google.com/file/d/18JPIjuPZPM2hnj5O18T4HeQnNil3fU9u/view?usp=sharing

https://drive.google.com/file/d/1RfIn2pIyr_zCyK18rQCX_J1Vr0OQ9Q4W/view?usp=sharing

https://drive.google.com/file/d/1l6M5xlKulDGQ2aFLz21TuNYY9s2LpFAd/view?usp=sharing

https://drive.google.com/file/d/1s86a0gtDE-15h0r3UORfMux3pT8Ns1ls/view?usp=sharing

According to your answer I tried to get the files which might be useful regarding annotation:

thanks a lot,  how can you do all of that within 24 hrs?  Congratulations!


text:
Reference variant detected as altered one in bam file
I received (from manufacturer) several .bam files and I used four callers (samtools, freebayes, haplotypecaller, deepvariant) to find some sequence variants. In obtained .vcf files, I took a closer look to some calls. I found interesting, homozygous one rs477033 (C/G Ref/Alt) with flag 'COMMON=0' and very low MAF. I also checked the .bam file (in IGV) and this rs itself in ensembl database. It looks like this:

The problem is: in this case ref and alt are switched. The G is reference (ancestral) allele and C is altered. Nevertheless in the .bam file it's interpreted in opposite way. Of course it's false positive variant - there's no mutation there at all. All variant callers returned it as high quality variant. My. bam files are full of them.
What went wrong in this case? Problem with mapping?
How can I fix this thing, to get only novel and rare true positive variants?

[link] https://drive.google.com/file/d/1eheG5MqS_a1xK5HSqsNJ0Tito_WdL_YV/view?usp=sharing

A reference sequence can obviously not represent all populations,  as shown above in the linked chapter (• page 30 •),  a reference allele may be the minor allele in the population you are studying ( in some cases )

The ancestral allele is not always the reference allele. Ref and alt aren't switched, ref is just rare in at least the population you're looking at (it appears to be rare in humans as a whole). So, the C is indeed reference and G is indeed altered, it's just that people homozygous for the G are the overwhelming majority. This is evidenced by the rs477033 variant being named NC_000001.11:g.165890393C>G. So, in short, nothing went wrong, all tools functioned correctly.
You'll notice a bunch of these over time, either make a database of them (or find one, someone may have made one) and filter them out or see if maybe the 1000 genomes reference genome has some of these modified (or maybe the T2T genome that was recently published).

I think that I found the solution and I'd like to add it here. It might be helpful for others. I just simply performed additional annotation of my .vcf files using 1KGenomes data. I've found it here: here Then each variant is described in INFO by AA which means ancestral allele. If this allele is different comparing to my REF I filter it out.


text:
Mixcr batch processing
I have a folder with 150 pairs of fastq files from an illumina sequence run. I need to precess them with mixcr, how can I do this in Bash with a single command?

What is your command for processing a single pair of FASTQ files?

You can use a for loop.

Here's an example of the command line that I'm using to run Mixcr for a pair of samples:  mixcr analyze amplicon --s hs --staring-material dna --5-end v-primers --3-end j-primers --adapters no-adapters --receptor-type trb --only-productive /Documents/Data/Input/S1_S1_L001_R1_001.fastq.gz /Documents/Data/Input/S1_S1_L001_R2_001.fastq.gz
The S1_S1_ component of the file name changes with each sample to S2_S2_  S3_S3_  all the way to 75.  I'm having a little problem setting the $variable for the file name in the loop.

One solution I thought is to try to iterate through the two types of files (R1 and R2) with a FOR loop, but I can't figure out how to do this.  I don't think a nested loop will do the trick.

Create a text file with two columns, one for R1 and the other for R2. Loop through this file. Shell loops are not a good thing to begin with, and nested loops are almost never required.

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

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

Thanks.  This is essentially the approach that I took.  My original problem actually was the result of trying to run a shell command that wasn't compatible with the Linux version of shell.  When I used $ batch instead of $ sh to initiate my script, all the recommendations worked. I should have started by providing the bit of information that I am working in Linux (Ubuntu 18.04 to be exact).

Ah - OK. Yes, you may want to use bash rather than sh if possible. Glad it's working for you.

I would also suggest using GNU parallel as described in our guides :
https://docs.milaboratories.com/mixcr/guides/generic-multiplex-bcr/#one-command-solution
Briefly:
ls /Documents/Data/Input/*_R1* | \
    parallel -j 4 \
    'mixcr analyze
    …
    {} \
    {=s:R1:R2:=} \
{=s:.*/:/Documents/Data/result/:;s:_R.*::=}'

That will take all fastq pairs from /Documents/Data/Input and put the result output files into /Documents/Data/result folder.


text:
PLINK linear association testing command with an interaction and 2 covariates term
I am trying to write the right command for my association testing using PLINK 1.9 but I could not find a clear answer so far and not sure how to wrtie it correct.
I have genotype variants called and in PLINK bed formatted files under geno_data folder.
I have phenotype.txt file. First column as "Sample_ID", second column as "Treatment" (1 for Placebo and 2 for CompoundX) and third column has the Response variable as "ResX", which has integer number between 0 to 9.
I also have PCs file as cov.txt. The first and second columns are the same and includes Sample_IDs. The column names are "FID" and "IID". The third, fouth and fifth columns have PCs,with the column names "PC1", "PC2" and "PC3".
The Plink documentation did not help me as I am new in this. Therefore I need some examples command to clarify for various models as described below.
I need the correct PLINK (1.9) command (if not possible PLINK2 is also OK) for the following linear models:
Model 1:
ResX ~ Variant * Treatment + PC1 + PC2 + PC3

Second command I need is for this model:
Model 2:
ResX ~ Variant + Treatment + PC1 + PC2 + PC3

Model 3:
ResX ~ Variant * Treatment

Note I had asked this question on Biostar 11 days ago but did not receive any answer or comment yet.
Thank you


text:
How to subset samples from a VCF file?
I have VCF files (SNPs & indels) for WGS on 100 samples, but I want to only use a specific subset of 10 of the samples. Is there a relatively easy way to pull out only the 10 samples, while still keeping all of the data for the entire genome?
I have a script that allows me to pull out regions of the whole genome for all 100 samples, so if I could do something similar but only put regions for the 10 samples that I want that would be ideal.

Welcome. Do you want random 10 samples or do you want to extract specific 10 samples? Also, do you talk about vcf files containing SNPs/indels ? What is the reason why you want to subsample your vcf file? We like details, details make questions understandable and answerable. You can [edit] your Q and add all the details.

Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the -s or -S option:

-s, --samples [^]LIST

Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. bcftools view is the exception where some
tags will be updated (unless the -I, --no-update option is used; see bcftools view documentation). To use updated tags for the subset in another command one can pipe from view into that command. For example:

-S, --samples-file FILE

File of sample names to include or exclude if prefixed with "^". One sample per line. See also the note above for the -s, --samples option. The command bcftools call accepts an optional second column indicating ploidy (0, 1 or 2) or sex (as defined by --ploidy, for example "F" or "M"), and can parse also PED files. If the second column is not present, the sex "F" is assumed. With bcftools call -C trio, PED file is expected. File formats examples:

sample1    1
sample2    2
sample3    2

or

sample1    M
sample2    F
sample3    F

or a .ped file (here is shown a minimum working example, the first column is ignored and the last indicates sex: 1=male, 2=female):

ignored daughterA fatherA motherA 2
ignored sonB fatherB motherB 1

Example usage:
bcftools view -s sample1,sample2 file.vcf > filtered.vcf
bcftools view -S sample_file.txt file.vcf > filtered.vcf

See the bcftools manpage for more information.

In addition to the answer from @gringer there is a bcftools plugin called split that can do this, but gives you the added ability to output single-sample VCFs by specifying a filename for each sample.
$ bcftools +split

About: Split VCF by sample, creating single-sample VCFs.

Usage: bcftools +split [Options]
Plugin options:
   -e, --exclude EXPR              exclude sites for which the expression is true (applied on the outputs)
   -i, --include EXPR              include only sites for which the expression is true (applied on the outputs)
   -k, --keep-tags LIST            list of tags to keep. By default all tags are preserved
   -o, --output DIR                write output to the directory DIR
   -O, --output-type b|u|z|v       b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]
   -r, --regions REGION            restrict to comma-separated list of regions
   -R, --regions-file FILE         restrict to regions listed in a file
   -S, --samples-file FILE         list of samples to keep with second (optional) column for basename of the new file
   -t, --targets REGION            similar to -r but streams rather than index-jumps
   -T, --targets-file FILE         similar to -R but streams rather than index-jumps
Examples:
   # Split a VCF file
   bcftools +split input.bcf -Ob -o dir

   # Exclude sites with missing or hom-ref genotypes
   bcftools +split input.bcf -Ob -o dir -i'GT="alt"'

   # Keep all INFO tags but only GT and PL in FORMAT
   bcftools +split input.bcf -Ob -o dir -k INFO,FMT/GT,PL

   # Keep all FORMAT tags but drop all INFO tags
   bcftools +split input.bcf -Ob -o dir -k FMT

So if you had the following samples file samples.tsv
sample1    sample1
sample2    sample2

You can run it and get the following
$ bcftools +split -S samples.tsv -o outdir in.vcf
$ ls
in.vcf sample1.vcf sample2.vcf

Without the second column, you would just get a single VCF with the two samples in it (as you would with view)

You can use the GATK's SelectVariants tool with the -sn flag.
E.g.

gatk SelectVariants -V input.vcf -R reference.fasta -sn Sample_01 -out sample.vcf

You may use the -sn flag several times so as to select several samples, or use it to point to a file containing a sample name on every line.

From the old VCF manual,
https://vcftools.sourceforge.net/man_latest.html
and
https://manpages.debian.org/testing/vcftools/vcf-subset.1.en.html
cat in.vcf | vcf-subset -r -t indels -e -c SAMPLE1 > out.vcf
comma-separated list of variant types ( indels ) could be obtained.

This is part of vctools, which is pretty much obsolete. I don't understand why you point to the ubuntu man page when vcftools has the closer-to-source manual.

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