text:
How do I include silent mutations in an oncoplot using maftools?
I have annotated files from annvar and I changed them into MAF by using
annovarToMaf from ‘maftools’. Then I used read.maf to read them into R and then I used merge_mafs to integrate them together inside the last file.
I had several data frames but it separated the synonymous mutations (maf.silent) from the non-synonymous mutations, and I want to create an oncoplot using the oncoplot function for both synonymous mutations from the non-synonymous mutations.
But once I use this command I don't know how to integrate both data frames inside of the MAF file.
How can I create the oncoplot using oncoplot for both synonymous mutations and non-synonymous mutations?

I have tried this, not work, the LOW impact SNP still missing
df.maf <- read.maf(maf = df, vc_nonSyn = c("Frame_Shift_Del","Missense_Mutation","Nonstop_Mutation","Silent","Splice_Region")) the reason I want to include silent SNPs is that I want to plot tumor mutation signature per sample. however, I got error due to the few number of SNPs for some samples.

This is late but I am putting it here for future reference. You can use the vc_nonsyn argument in the read.maf function to manually list the variant classifications you want to plot later. It will look like this:
df.maf <- read.maf(maf = df, vc_nonSyn = c("Frame_Shift_Del","Missense_Mutation","Nonstop_Mutation","Silent","Splice_Region"))


text:
Allele Count and Allele Frequency in VCF files
I'm working in bioinformatics, but my computational skills far outstrip my knowledge of biology or genomics. So forgive the noobish question.
According to the VCF specification, the INFO column for each coordinate [CHROM, POS] can contain values for
AC : allele count in genotypes, for each ALT allele, in the same order as listed
AF : allele frequency for each ALT allele in the same order as listed: use this when estimated from primary data, not called genotypes
AN : total number of alleles in called genotypes
For example, I'm working with the the VCF file you can download from Kaviar.
I think I understand how to calculate those numbers from raw sequences from a pileup. If for a given coordinate [CHROM, POS], $N_0$ is the number of hits for the reference allele, $\textrm{AC}_i \equiv N_i$ is the number of hits for the $i^{th}$ variant (maybe only one), and AN is the total number of alleles counted: $\textrm{AN} \equiv \sum_0^n N_i$, then the allele frequency for the $i^{th}$ variant is
$$\textrm{AF}_i = \frac{\textrm{AC}_i}{\textrm{AN}} = \frac{N_i}{\sum_0^n N_i}$$
So the phrase "total number of alleles in the dataset" would mean "total number of hits to coordinate [CHROM, POS] for any allele. And by "hit" I mean an aligned sequence fragment that includes that coordinate.
But I may not be taking into account the caveat for AF: "use this when estimated from primary data, not called genotypes."
Am I correctly understanding AF, AC, and AN in VCF files?

No, that is not quite right. specifically, this part:

And by "hit" I mean an aligned sequence fragment that includes that coordinate.

To get data on how many reads contributed to calling a specific allele at a given position, one would use the DP field, for example, which isn't included in the Kaviar VCF.
The AC, AN, and AF fields in a VCF file are meant to be used to show the allele's frequency in the context of all of the samples used when making that VCF file.
Investigation process
This is one line from the Kaviar vcf file:
1       10230   rs200279319     AC      AA,A    .       .       AF=0.0002654,0.0048525;AC=7,128;AN=26378;END=10231

If we look at the README attached with the Kaviar vcf file, it says,

"Kaviar is a compilation of SNVs, indels, and comlex variants observed in humans, designed to facilitate testing for the novelty and frequency of observed variants."

How many humans though? If we go to the Kaviar website it says:

Kaviar contains 162 million SNV sites (including 25M not in dbSNP) and incorporates data from 35 projects encompassing 77,781 individuals (13.2K whole genome, 64.6K exome).

If we look at the Kaviar VCF header we see that the fields are defined:
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele Count">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in data sources">

If we look at these fields again, AF=0.0002654,0.0048525;AC=7,128;AN=26378, the AN=26378 tells us that for this locus, the Kaviar dataset contains information from 26378/2=13189 humans. That's the 13.2K whole genomes the Kaviar website was talking about. This site I'm looking at must be a non-coding region, because there is no exome data (protein-coding regions).
For this same site, AC=7,128 means that out of 13189 humans (26378 haplotypes), we saw the variants 7 and 128 times. Out of 26378 haplotypes, this is a frequency of AF=0.0002654,0.0048525.

Thank you so much for that clear and detailed answer!

source: Stack Exchange

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

sort: votes


text:
Fitting branch models with godon
Is it possible to fit branch models (Yang and Nielsen, 1998) in godon? That is, the model which averages over all sites, as opposed to the branch site models?

For those VTCing this - Bioinformatics.SE is the specified help forum for the godon package. If you've got a problem, take it up with the meta

The votes to close are because your question lacks a lot of detail.

And yet, it has all of the information necessary for the question to be answered.

Voted to open, because it's exactly the type of questions we aim for - concise, concrete and answerable.

Branch-models are currently not implemented in godon. In the PAML manual Ziheng Yang writes:

This model is very parameter-rich and its use is discouraged.

This was one of the reasons I did not implement it.
Disclaimer: I am the author of godon.

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


text:
How to extract all peptide sequences for a UniProt protein entry with an API?
I know how to generally pull down information for a UniProt entry using the REST API, for example:
https://www.uniprot.org/uniprot/?query=id:Q9UJL9&sort=score&columns=id,go(molecular%20function)&format=tab
I'm looking at this entry as an example. I'm able to pull down most of the attributes I want with a variation of the API query above, however, I'm stuck on one part.
For all human proteins, I want to pull down the peptide sequences that come from that protein.
In the example above, it would be the sequences of the peptides associated with this section:

Can this be done with the API or do I have to come with a scraping method for all human proteins where I (1) get the protein sequence of each entry, (2) get this table, (3) calculate where the peptides start and end from the sequence in step 1?
I feel like this alternative will be much slower, so I want to check I'm not missing anything obvious.

Downloading the whole of Uniprot and filtering for humans is probably a nicer solution. It is in XML format, which is a bit of a pain, but is a goldmine of data, which you may need —totally depends on your plans and how they may evolve. E.g. I have written my own Python XML Uniprot parser due to my needs and I have a human scraped dataset that is heavily focused on my project which expanded with time.

If downloading the complete database is not an option, you may able to use the instructions from this help page (based on the gff format for any given query):
https://www.uniprot.org/help/download_subsequences


text:
Finding SNP/Indel Markers of Interspecific Variation

This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.
Goal:
Find SNP and indel alleles which are fixed between but lack polymorphism within two closely related species (I'm looking for the BLUE region of the Venn diagram)
Desired Outputs:

All SNPs/Indels of Species A which are exclusive to and lack variation within Species A

All SNPs/Indels of Species B which are exclusive to and lack variation within Species B

Scratch Materials:

Species A reference genome

Species B reference genome

Illumina short reads for Species A

Illumina short reads for Species B

My current plan is to...

Map Species B to Species A:

minimap2 -ax asm5 referenceA.fasta referenceB.fasta > BtoA.sam

Find all SNP/Indel alleles of Species B relative to Species A (Represents the entire circle in the above Venn diagram):

bcftools mpileup -f referenceA.fasta -s BtoA.sam | bcftools call -o B_compared_to_A.vcf

Find all SNPs/Indel of Species B reads relative to Species B reference (Represents the red region of the Venn diagram):

bcftools mpileup -f referenceB.fasta -b bam_listB | bcftools call -o B_compared_to_B.vcf

Perform a liftover on "B_compared_to_B.vcf" to correct for coordinate differences between Species A and Species B references (I haven't completely figured out how to do this - I want to make sure I have the general pipeline nailed down before I get into the weeds)
Write a program which extracts all loci found in "B_compared_to_A.vcf" but NOT in "B_compared_to_B_liftedOver.vcf" (This should yield the BLUE region of the Venn Diagram
Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

That seems like a valid approach to me.

I might be missing something, but does this approach not require that the coordinate systems of referenceA.fasta and referenceB.fasta be the same? I would expect a liftover step to map variant coordinates between the two. For example, the BAM headers in bam_listA will not match the fasta index of referenceB.fasta. Maybe these tools don't care about that, I don't know that off the top of my head. but I worry about coordinate frameshifting.

At the risk of  sounding dumb, is there a reason that you don't use reciprocal whole-genome alignment via e.g. minimap2 to identify the interspecific differences, which you can then check with reads in each species separately? i don't have a specific workflow fleshed out, but i would hope that the reference assemblies were good enough to get the (vast) majority of things that would be picked up by reciprocal shotgun read mapping.

@MaximilianPress could you link me to an article describing coordinate frameshifting? I am relatively new to bioinformatics and haven't encountered this concept before. I will also give minimap2 a look, I haven't heard of this software before.

not really an article-level concept (and others would use different term probably), it's just that the two genomes probably don't have the exact same coordinates (e.g. chr1 is 10.5 Mbp in one genome, and 10.6 Mbp in the other genome). If this is true, then when you call a variant in one genome, those coordinates will usually differ from the coordinates for the homologous variant in the other genome. Stated differently: it is very unlikely that two species' genomes can be aligned without gaps. Wherever there is a gap, variant coordinates will shift.

@MaximilianPress I apologize for dredging up a nearly month-old conversation, but I am still struggling with this problem. When you mentioned "reciprocal whole-genome mapping" with minimap2, what did you intend for me to do with the reads once I have mapped the Species A reference genome to the Species B reference genome? I understand you didn't have a pipeline worked out, but I'm trying to figure out the logic of the approach.

I will post as a short answer, easier to format etc.

Ok, with invaluable aid from @Maximillian Press, I've figured out how to find alleles which are variable between two species but fixed within one species.
Here is how you find SNP/Indel loci which are distinct between species A (hereafter sA) and species B (hereafter sB) but are fixed WITHIN sA:
# Maximillian Press's Contribution: Map sA reference to sB reference using minimap2, then convert to bam.

minimap2 -a sB_ref.fasta sA_ref.fasta > AtoB.sam
samtools view -bh AtoB.sam > AtoB.bam

# Sort bam by coordinates, then use bcftools to generate a pileup and call all alleles of sA relative to sB

java -jar $EBROOTPICARD/picard.jar SortSam I=AtoB.bam O=sorted_AtoB.bam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR
bcftools mpileup -O b -o AtoB.bcf -f sB_ref.fasta sorted_AtoB.bam
bcftools call -m -v -o AtoB.vcf AtoB.bcf

 # Call all alleles of sA relative to sA.
# This represents the SNP/Indel loci which are variable WITHIN sA and thus ought to be discarded from the final output.
# "bam_listA.txt" contains names of all the coordinate-sorted bam files for species A
bcftools mpileup -O b -o AtoA.bcf -f sA_ref.fasta -b bam_listA.txt
bcftools call -m -v -o AtoA.vcf AtoA.bcf

I then use the "sed" function to replace chromosome names with their analogs (happy to give more details if anyone else needs to do this) and also removed the VCF headers prior to making a C++ program which extracts all the rows found in AtoB.vcf but not in AtoA.vcf. My final output was thus all SNP/Indel loci which vary between sA and sB but NOT within sA.
Happy to share the C++ program if anyone else needs to do this.

Problem statement
You want to find a set of SNPs in populations of two species, for each of which you have a distinct reference genome. You then want to know which of those SNPs are shared between the two genomes (there is variation in each species), and which ones are private to only one of the species (in the other genome there is no variant at that location).
To accomplish this task, you need to first find SNP variation in each genome, and then figure out which SNPs are the same between the two genomes. In other words, line up two VCF files for different genomes and figure out which SNPs are the same between them. In comparing genome assemblies, this is sometimes called a liftover operation.
Genomes can have different coordinate systems.
Let's first think about what exactly you are looking at with 2 different reference genomes. For example, here are outputs of samtools faidx for 2 different S. cerevisiae strains, S288C and Kyokai-7:
# S288C
ref|NC_001133|  230203  95  60  61
CHR2.19980913   813142  234173  60  61
CHR3.19980521   315339  1060907 60  61
CHR4.19990210   1531929 1381540 60  61
CHR5.19970727   576870  2939039 60  61
CHR6.19970727   270148  3525562 60  61
CHR7.19970703   1090936 3800252 60  61
CHR8.19970727   562638  4909411 60  61
ref|NC_001141|  439885  5481523 60  61
CHR10.19970727  745440  5928778 60  61
CHR11.19970727  666448  6686681 60  61
CHR12.19970730  1078172 7364277 60  61
ref|NC_001145|  924430  8460517 60  61
ref|NC_001146|  784328  9400452 60  61
CHR15.19970811  1091283 10197892    60  61
ref|NC_001148|  948061  11307461    60  61
ref|NC_001224|  85779   12271444    60  61

# Kyokai first 10 lines
gi|347729946|dbj|BABQ01000001.1|    5491    137 80  81
gi|347729945|dbj|BABQ01000002.1|    4704    5834    80  81
gi|347729944|dbj|BABQ01000003.1|    203046  10734   80  81
gi|347729943|dbj|BABQ01000004.1|    1036    216456  80  81
gi|347729942|dbj|BABQ01000005.1|    126714  217642  80  81
gi|347729941|dbj|BABQ01000006.1|    121900  346077  80  81
gi|347729940|dbj|BABQ01000007.1|    20385   469638  80  81
gi|347729939|dbj|BABQ01000008.1|    243458  490415  80  81
gi|347729938|dbj|BABQ01000009.1|    189876  737054  80  81
gi|347729937|dbj|BABQ01000010.1|    81973   929441  80  81
# ...
# goes on for many thousands of lines

What you will notice is that there are significant differences between the 2 assemblies. The first is almost all in chromosomes, the second is largely mapped to various chromosomes but exists in thousands of tiny contigs. There is no way to just map reads to them, call variants, and assume that you (or bcftools) can figure out what is going on.
This is possibly an extreme example, but even in the case of two genomes assembled to chromosome pseudomolecules, there is basically no chance that those two chromosomes are on the same coordinate system or have similar metadata. That means again that bcftools probably won't be able to figure it out, because a SNP at position 112929 in genome A will be homologous to position 110201 in genome B.
Liftover: a way to solve the mapping problem between 2 genomes.
You can create a mapping between the two genomes. This is 100% a solved problem in bioinformatics, and it is more or less the same as doing liftovers between different assembly versions of the same genome. I think that you can more or less just use the UCSC LiftOver tool to convert SNP coordinates, unless I am missing something. That how-to guide that I linked is going to give you much better instructions than I can give you, here is a more minimal version that might be simpler.
The only extra suggestion that I have is that the guide uses BLAT to align the two assemblies, and minimap2 is a much faster tool. You may need to fight the input formatting a little bit, but if you run into runtime issues you might make that change. This is why I suggested it in the comment, though I wasn't aware that the liftover docs were so complete.
Short pipeline:

Align reads (individuals of species A to reference A, individuals of species B to reference B).
call variants-->get VCF files using preferred workflow for each of A and B independently.
compute liftover mapping (chain file, usually) from reference A to reference B.
use computed liftover mapping to actually lift VCF over from one genome to another.
align A and B references, call SNPs between them (possibly can be done from (3) output). Ensure that you are using as reference the liftover target genome so that the VCF is informative.
For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.

Thank you for your detailed response! I still am not sure at which point in the pipeline I perform variant calling. Do I use the liftover output in lieu of the Species A and Species B reference genomes in the SNP-calling steps to get the SNPs of each species? Here's an example of what I'm thinking of:
bcftools mpileup -f liftover  -b bam_listA | bcftools call -o ACompatibleCoordinates.vcf bcftools mpileup -f liftover -b bam_listB | bcftools call -o BCompatibleCoordinates.vcf
I would then search the output VCF files for alleles which are unique to each species

I have added a short suggested workflow that indicates what I saw as the best option. I am not sure that I totally understand your suggestion, it looks like you are passing a liftover to bcftools but I don't totally understand what that is supposed to look like? You could pass a FASTA which is the output of a liftover to bcftools I guess, or just use it as a reference in alignment in the first place. Those would likely also work, though I don't have a great sense for which way is best.

Ok, thank you! I see a problem with the proposed workflow. In my understanding, to be called as an SNP/Indel, the locus in question must have different alleles in the reference and sam/bam. I am trying to find all SNPs/indels which are exclusive to Species A and all SNPs/indels which are exclusive to Species B (i.e. completely lack variation within species). If I call SNPs/Indels of A relative to A, won't the resulting VCF file not include alleles which are 100% fixed in A? I am about to post an answer to my own question which solves this issue, though I would appreciate critique

@annabelperry in this case, can you not align the two genomes to call such variants, and then confirm post hoc that they do not have variation in the BAM files? e.g. minimap2 -a genomea.fa genomeb.fa | samtools view -bh | bcftools etc., and then scan those against whichever VCF is the reference? that should be quite simple.

Ah, thank you @Maximilian Press! I will try that and upvote your answer if it works. Else, I will be back here shortly haha


text:
How can I retreive specific qualifiers from genbank files in R?
I'm trying to extract the host field from the features table using the biofiles package. From what I understand, the function filter (from the same package) could be of help but I can't extract the information of interest. This is the simplified series of steps that I'm following:
library(biofiles)
x_search <- reutils::esearch(db = "nuccore", #retmax = 9999,
             term = "((bacteria[ORGN]) AND host[ALL]) AND Loripes[ALL]")

x_fetch <- reutils::efetch(x_search[1:3], "nuccore",rettype = "gbwithparts",
            retmode = "text")

x_rec <- biofiles::gbRecord(x_fetch)

x_ft <- biofiles::getFeatures(x_rec)

x_ft[1]

The last output should show the list of features for that entry, including the source with the qualifier "host", which I would like to retrieve for all the ids.
EDIT: I'm open to solutions outside the biofiles package.

Easy with biopython if Ur open to answers outside of R

@LiamMcIntyre Sure! I would like to know your suggestion.


text:
Filtering samples by FORMAT attributes using bcftools
The -e and -i options of the bcftools filter command appear, by default, to only allow for including or excluding sites.  E.g., -e 'FMT/DP < 10' removes sites where any sample has DP < 10, and -e 'MEAN(FMT/DP) < 10' removes sites where average depth across samples is < 10.
I would like to perform effectively similar filtering commands, but in a way that includes or excludes samples, instead of sites.  E.g., exclude all samples that have DP < 10 at any site, or exclude all samples with an average depth across all sites < 10, or exclude all samples that have a genotype quality of < 20 at 10% of genotyped sites.
Such commands can very easily be written to apply to sites, but there does not seem to be a straightforward way of applying such filtering expressions to samples in bcftools.
One option mentioned in the manual is to use the double || operator, as described in the man page.  As described in my comment on a Github issue, however, this does not seem to work as a filtering operation for including or excluding entire columns in a VCF file, but rather only for selecting more or less individual records (e.g., as illustrated in the example on this tutorial page).  In addition, it would seem limited in functionality (e.g., functions such as MEAN(FMT/DP) < 10 would still be evaluated across samples, not sites, when combined with a || operator).
The only other method mentioned in the manual is to use the SMPL or s prefix to functions which evaluate FORMAT fields in expressions passed to filter.  This also does not appear to work as I would expect.  I'm not sure I understand what it is doing, but it definitely isn't removing columns from the VCF file.  Instead it seems to be using as a metric calculated across sites on a per sample basis as a filtering criteria applied to sites (but I could just be confused):
nSites_unfiltered=$(bcftools view -H unfiltered.vcf.gz | wc -l)
nSamples_unfiltered=$(bcftools query unfiltered.bcf.gz -l | wc -l)

nSites_DPfilter=$(bcftools filter unfiltered.bcf.gz -e 'sMEAN(FORMAT/DP)<5' | bcftools view -H | wc -l)
nSamples_DPfilter=$(bcftools filter unfiltered.bcf.gz -e 'sMEAN(FORMAT/DP)<5' | bcftools query -l | wc -l)

echo "scale=5; $nSites_DPfilter / $nSites_unfiltered" | bc #0.17349 - sites are removed, even though the filter is being applied to samples
echo "scale=5; $nSamples_DPfilter / $nSamples_unfiltered" | bc #1.00000 - no samples removed, despite, again, the filter being applied to samples

nSites_normalDPFilter=$(bcftools filter unfiltered.bcf.gz -e 'MEAN(FMT/DP) < 5' | bcftools view -H | wc -l)
echo "scale=5; $nSites_normalDPFilter / $nSites_unfiltered" | bc #0.83121 - fewer sites removed than the sMEAN filter

The .bcf.gz passed to these commands is here: https://www.dropbox.com/s/fwj1snabsb763sn/unfiltered.bcf.gz?dl=0
Its .csi index file is here: https://www.dropbox.com/s/jzdq6y91nkencmy/unfiltered.bcf.gz.csi?dl=0
In any event, the only option I have been able to come up with that actually excludes or includes entire sample IDs is to extract summary statistics using bcftools stats, export the desired statistics from the PSC, Per-sample counts table, and perform my own exploratory analysis outside of bcftools.  Then I can use bcftools view -S to explicitly remove samples by their ID, which I have identified manually as failing some thresholding criteria.  This also is obviously cumbersome, and only works for FMT fields that are included in stats (e.g., not GQ or PL).
I am not sure if I am missing some basic functionality of bcftools, or if this is more of a feature request.  Any information would be appreciated.
Thanks!

You may want to take a look at vcftools' genotype filtering options or the custom JEXL expressions that GATK's SelectVariants allows. Are you sure your sample-filtering operations won't produce a variable number of columns for each row? Also, remember that bcftools filter will not remove samples no matter what you try. It's a utility to filter sites, and it can only filter sites.

Thanks very much.  Since every site has genotypes for every sample (if you include missing genotype entries), I'm not sure I understand how removing certain columns in their entirety would create any problems in terms of uneven dimensions.  Indeed, this is the function of bcftools view -s ^<list.of.sample.IDs>, if I am not mistaken.  vcftools also only allows for filtering samples based on ID.  Thank you for the reference to SelectVariants.  Do you have a link to a manual, or example code to address any of the example use-cases referenced in my question?  Thanks!

Not really, no. I have rarely ever used SelectVariants when it wasn't already implemented in a pipeline (GATK tools can be a pain to use at times). If push comes to shove, it might be worth converting the VCF to MAF format so the data is at individual level.


text:
Kallisto psudoaligner quant not giving me expected output
I am trying to run this kallisto (latest) quant and all I get are abundance and json files.
kallisto quant --plaintext -t 4  -o rf-stranded.nob.1K1_L7.0 -i mm9_ref.idx  --rf-stranded 1K1_L7.0 2K1_L7.0

the results files in stderr is:
[quant] fragment length distribution will be estimated from the data
[index] k-mer length: 15
[index] number of targets: 23
[index] number of k-mers: 363,704,628
[index] number of equivalence classes: 1,945,936
[quant] running in paired-end mode
[quant] will process pair 1: 1K1_L1.0
                             2K1_L1.0
[quant] finding pseudoalignments for the reads ... done
[quant] processed 22,897,424 reads, 4,476,549 reads pseudoaligned
[quant] estimated average fragment length: 24.2427
[   em] quantifying the abundances ... done
[   em] the Expectation-Maximization algorithm ran for 52 rounds

I get an abundance file with:
target_id       length  eff_length      est_counts      tpm
chr1    197195432       1.89745e+08     253015  21.3249
chr10   129993255       1.31621e+08     210320  25.5544
chr11   121843856       1.38776e+08     263606  30.3773
chr12   121257530       1.23195e+08     143568  18.6371
chr13   120284312       1.22484e+08     132969  17.3613
chr14   125194864       1.20926e+08     112711  14.9058
chr15   103494974       1.07866e+08     109335  16.2102
chr16   98319150        9.70075e+07     109045  17.9767
chr17   95272651        1.06375e+08     169166  25.4322
chr18   90772031        9.23542e+07     121990  21.1241
chr19   61342430        7.23377e+07     157142  34.7407
chr2    181748087       1.81791e+08     244870  21.5414
chr3    159599783       1.45978e+08     209257  22.9247
chr4    155630120       1.61978e+08     226893  22.4014
chr5    152537259       1.65672e+08     194951  18.8186
chr6    149517037       1.43243e+08     246652  27.5373
chr7    152524553       1.82795e+08     197527  17.2812
chr8    131738871       1.50556e+08     203954  21.6643
chr9    124076172       1.30446e+08     222639  27.295
chrX    166650296       1.40278e+08     176259  20.0942
chrY    15902555        5.76797e+07     123.462 0.0342311
mtDNA   16299   11804.8 737285  998819
rDNA    249368  720854  33272.5 738.159

run_info.json:
        "n_targets": 23,
        "n_bootstraps": 100,
        "n_processed": 22897424,
        "n_pseudoaligned": 4476549,
        "n_unique": 3641646,
        "p_pseudoaligned": 19.6,
        "p_unique": 15.9,
        "kallisto_version": "0.46.1",
        "index_version": 10,
        "start_time": "",
        "call": "kallisto quant --plaintext -t 4 -b 100 -o rf-stranded.1K1_L1.0 -i mm9_ref.idx --bias --rf-stranded --gtf UCSC_genes.gtf --pseudobam --genomebam 1K1_L1.0.gz 2K1_L1.0.gz"
}

When I ran it with boot strapping (-b), -gtf and --pseudobam --genomebam it gave me a bam  file.
44879990 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
8078785 + 0 mapped (18.00% : N/A)
44879990 + 0 paired in sequencing
22439995 + 0 read1
22439995 + 0 read2
8077160 + 0 properly paired (18.00% : N/A)
8077160 + 0 with itself and mate mapped
1625 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

I'm not sure what is going wrong, can someone help?
NOTE I use a 15nt kmer, and this about 1/2 of what is suggested (31 I think), is that a problem?

Thanks I think I know what is going on.  I probably should give a file of transcripts to the index step and not a chromosome file.  This is why it is giving my only chromosome abundances.

That makes sense. I would suggest posting that as an answer to your own question and then accepting it, to help others.

OK but I first will try to verify.

I probably should give a file of transcripts to the index step and not a chromosome file
Not probably, you must. Kallisto strictly pseudoaligns against a transcriptome, that is its mode of operation.

As you surmised, kallisto requires a transcriptome fasta file, not a genome fasta file.


text:
Genome QC + Assembly Pipeline semantics
I’m trying to create a pipeline for genome assembly. How best can I “redirect/pipe” from existing fasta files (or files in general) to other steps of the pipeline?
I was thinking of going from the SRA Download step (prefetch and fasterq-dump) to fastqc, then to trimmomatic or bbduk to remove the adapters, one more fastqc check and then to SPAdes for the assembly.
In my (amateur) mind it would be cool to be like:
“SRRXXXXXX” | prefetch | fasterq-dump | fastqc | bbduk.sh | fastqc (second) | spades.py
Any and all help is greatly appreciate! I know this stuff is pretty basic but I am so new to everything it’s painfully obvious to even me.

Maybe take a look at workflow management systems, like

snakemake
Nextflow
Toil
Cromwell

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

I understand the desire to pipe everything like that but I don't think it will work, even for such a small workflow. And as your pipe grows it'll make it harder to troubleshoot/debug.
These days everyone uses workflow languages to write pipelines. It's a bit more work to get set up but worth it in the long run. I use WDL. Snakemake is very popular. BASH will of course work for small pipes but not with | for every step.

Directly piping will definitely not work, as Liam McIntyre suggests. However, if you don't feel up to learning a workflow language just to do this, it is perfectly possible to just put the various commands in a bash driver script (or whatever you're comfortable with).
Each of those tools will take parameters to define output filenames, or have systematic ways of writing FASTA files and defining inputs. You will have to look at the command line interface of each of those tools to figure it out (for example, sratools here, bbduk.sh here and spades here). an example would be as follows:
#!/usr/bin/bash
SRA=$1

# get reads
prefetch $SRA  # writes to a known path locally
fasterq-dump $SRA --split-3  # looks for path_to_sra_file and outputs it. assumes you are splitting.

# trim reads
bbduk.sh in1=$SRA''_1.fastq in2=$SRA''_2.fastq out1=$SRA''_clean_1.fq out2=$SRA''_clean_2.fq

# re-QC the files- this step does not feed into the other steps!!
fastqc $SRA''*.fastq $SRA''*.fq  # runs fastqc on all fastq files

# run assembly
spades.py -o my_assembly --pe1 $SRA''_clean_1.fq --pe2 $SRA''_clean_2.fq

and so on.
I suppose that from a formal perspective it's better to use a workflow engine as Liam McIntyre says, but with such a straightforward workflow I'd say start simple and get more elaborate only if you need to; you'd spend 10x as long just figuring out the workflow language. This is ok if it's a goal to learn it or if you already know it, but if you just want some results putting it into a script is fine.

The other virtue of putting everything in a script is every command line option is documented.

All workflow management systems I have ever worked with also document all parameters (and much more) meticulously.


text:
Does anyone know how I can convert DNA code into FASTA for this TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC
TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC this is the DNA code Sequence. I have tried converting each into the amino acid and ended up with this LKHWMNEKPCFATPQX but I was told that this isn't right.

So you googled what fasta format looks like for DNA, and you concluded from what you read that it meant to translate the DNA into protein?

FASTA is a sequence format, and it is used for any sequence, including DNA and protein sequence. Translating DNA to protein sequence does not require FASTA or any format, for that matter - the sequence itself is enough. Please clarify your question.

This is your nuc seq as a fasta:
>sequence_1
TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC

To translate it use an online tool (ie https://web.expasy.org/translate/) or do it in biopython:
from Bio import SeqIO
for record in SeqIO.parse("/path/seq1.fa", "fasta"):
    print ('nuc', record, sep='\n')
    print ('prot', record.translate(), sep='\n')

Which gives:
nuc
ID: sequence_1
Name: sequence_1
Description: sequence_1
Number of features: 0
Seq('TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC', SingleLetterAlphabet())
prot
ID: <unknown id>
Name: <unknown name>
Description: <unknown description>
Number of features: 0
Seq('LKHWMNEKPCFATPQ', ExtendedIUPACProtein())

Note, this defaults to the human codon table. Other codon tables are available with the parameter 'table'.


text:
Diploid Consensus Sequences
My goal is to generate a diploid consensus genome. By "diploid consensus", I mean that I want to merge all supporting reads for a genomic region into a single haplome such that the resulting sam file contains TWO sets of reads, one from each haplome.
I have Illumina short-read files from the NCBI's SRA. These short-reads were obtained from a diploid species.
I am worried that, by simply using Picard's MarkDuplicates with the "REMOVE_DUPLICATES" option set to "TRUE", I will retain only the homolog with the highest quality and output a haploid version of my genome.
Here is my current workflow:

## Map Illumina reads to the reference genome. Include read IDs
bwa mem -M -R '@RG\tID:BioSample.LibraryName\tSM:BioSample\tLB:LibraryName\tPU:SingleFlowcell.SingleLane\tPL:ILLUMINA\t' reference.fasta SRR[number]_1.fastq SRR[number]_2.fastq > SRR[number].sam

## Sort Illumina reads by coordinate
java -jar $EBROOTPICARD/picard.jar SortSam I=SRR[number].sam O=sorted_SRR[number].sam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR

## Remove duplicate reads from sam file
## *Note: How do I remove sampling error reads?
java -jar $EBROOTPICARD/picard.jar MarkDuplicates REMOVE_DUPLICATES=true REMOVE_SEQUENCING_DUPLICATES=true I= 'sorted_SRR[number].sam' O= nodup_SRR[number].sam METRICS_FILE= metrics_SRR[number].txt TMP_DIR=$TMPDIR

## Convert duplicate-free sam file to bam file
samtools view -S -b nodup_SRR[number].sam > nodup_SRR[number].bam

## Re-align reads around indels to reduce slippage.
gatk LeftAlignIndels -R reference.fasta -I nodup_SRR[number].bam -O aligned_SRR[number].bam

Does anyone know of a program which will whittle my sam files down to ONLY two homologs per region?

What you are asking for doesn't make sense. How on earth is any software supposed to determine if a discrepancy is a polymorphism or an error if you remove all the supporting reads?

@swbarnes2 It's not that I want to simply remove the supporting reads - I should clarify my intent. My end goal is to use the consensus genome as input in a homebrew microsat-calling software. I don't want this software to call sequencing errors as microsats, so I want to merge all reads of the same sequence into a single consensus read based on base quality scores, generating two sets of sequences which are the pipeline's best guess at the individual's haplomes. I will edit my original post to make this clear.

@swbarnes2 I want to do what is described in this article, but I don't have PacBio long-reads to compliment my Illumina short-reads: https://doi.org/10.1093/bioinformatics/bty279

FWIW, I do not think that removing lines from a bam file is what you want to do at all.  I think you need to be calling variants with phasing data, and using that vcf to make a modified fasta.  A quick google turned this up  https://www.biostars.org/p/225262/

This is a little hard to answer as is, could you give an example of the kind of data you want to get out? Do you want to phase the haplotypes, like in the linked paper? I think swbarnes2 has the right of it, that calling variants is the way to go. Also- duplicates are usually defined as reads with the same coordinates; unless you have very high coverage this is unlikely to remove more than a very small fraction of informative reads (they are more likely PCR or optical duplicates, i.e. artifacts).

@MaximilianPress Ah I have failed to clarify my end goal. I am designing a software for microsatellite detection from short-reads (no such software currently exists which meets my needs). I only want to call variants in "real" reads rather than duplicates. However, I have found an alternative to diploid consensus creation - I am instead going to run my software on short-reads (including duplicates) and use variant frequencies to filter out false positives.

How to separate reads, one for each reference haplotype.

Map all the reads against the reference genome
Use GATK to call SNPs
Use GATK tool [FastaAlternateReferenceMaker][1] to build the alternative reference
You now have the two "haploid references" (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)
Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called [BayesASE][2]

In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved.
[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker
[2]: https://github.com/McIntyre-Lab/BayesASE


text:
How to filter out partially similar strings from two lists and make one list?
I have following two lists.
List1
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100008029/0_42112 id=244006
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100008044/0_60161 id=24400
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199
m64041_200717_231916/100008373/4102_50601 id=2440234

List2
m64041_200717_231916/100664277/0_42364
m64041_200717_231916/100665802/0_46854
m64041_200717_231916/100928073/0_11804
m64041_200717_231916/10094492/0_44157

These are sequence read ids from a genome file. List1 contains all the ids from the List2 but with some id numbers (id=2440). So the sequence ids are only similar to some extent. I want to compare these two files and I need to filter out List2 sequence ids from List one (with id numbers) as below.
Output.
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199

I will be grateful if someone can help me with this issue.

It sounds easy but I don't really understand exactly what you want to compare. Can U expand?

For another example imagine there are two lists. List1.
UnitedStates id=12,
UnitedKingdom id=15,
Spain id=20,
France id=25,
List2
UnitedStates,
Spain
I want to filter out UnitedStates id=12 and Spain id=20 using the List2. List1 and List2 are two files.
Output should be,
UnitedStates id=12,
Spain id=20

ah, so only keep whats in list 2?

Yeah. list2 will be the query.

You really need to learn some python :)
with open('/path/list2.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip())

with open('/path/list1.txt', 'r') as fin:
    with open('/path/list3.txt', 'w') as fout:
        for line in fin:
            if line.strip().split(' ')[0] in to_keep:
                fout.write(line)

Thank you very much. It worked. Yes, I need to learn python.

And another question. Do you know a way to filter out/remove those sequences from the original fasta file using the id list created from the above code? I have a piece of code I used earlier to remove sequences using a list. But since these ids contains a space, It's not working now. I tried to change it, but it didn't work. It'll be an easy thing for you,but I couldn't figure it out. Can you please help?

Sure. Post a new question with an example plz

Why would you hardcode two file names instead of just reading argv for the input file and printing to stdout so you can then redirect?

He asked for files

There is no need for a highlevel language in such a simple case, an awk oneliner does just fine:
$ awk -F " " 'FNR==NR {a[$1]; next} $1 in a' List2 List1
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199

This can be done in python easily. If you have different text files then you can load them into a python list and follow from step 2
# Step 1
list1 = ["m64041_200717_231916/100664277/0_42364 id=2440",
"m64041_200717_231916/100008029/0_42112 id=244006",
"m64041_200717_231916/100665802/0_46854 id=2440067",
"m64041_200717_231916/100008044/0_60161 id=24400",
"m64041_200717_231916/100928073/0_11804 id=2440152",
"m64041_200717_231916/10094492/0_44157 id=2440199",
"m64041_200717_231916/100008373/4102_50601 id=2440234"]

list2 = ["m64041_200717_231916/100664277/0_42364",
"m64041_200717_231916/100665802/0_46854",
"m64041_200717_231916/100928073/0_11804",
"m64041_200717_231916/10094492/0_44157"]

# Step 2
split_list1 = [ele.split(" ")[0] for ele in list1 ]
mapdict = {split_list1[i]:list1[i] for i in range(len(list1))}
common = [value for value in split_list1 if value in list2]
required = [mapdict[x] for x in common]
print(required)


text:
Books on application of Machine Learning in Bioinformatics
I'm looking for books on applications of machine learning and statistical data mining in bioinformatics with example codes provided in R and/or python.
I came across the following two books: the first being a quite good one, provides the code in MATLAB, which I'm not very familiar with. The second one too is a good book with examples in R.

Feature Selection and Ensemble Methods for Bioinformatics by Oleg Okun

Data Science and Predictive Analytics: Biomedical and Health Applications using R by Ivo Dinov

Can people suggest other books based on their first hand experiences?
Thanks

There are plenty of books when we talk about bioinformatics but only a few of them are available for application of machine learning in bioinformatics. From these selected books I would suggest that you go ahead with Statistical Modelling and Machine Learning Principles for Bioinformatics Techniques, Tools, and Applications. Its a recent book published in 2020 and covers most of the important topics in bioinformatics along with their ML applications. Apart from that you can also refer Machine Learning in Bioinformatics it’s quite old but still it covers most of the concepts with in-depth explanation. Along with that you can refer to this amazing book Computational Genomics with R available free online. This will cover both theory and code for different ML algorithms (both supervised and unsupervised) with their application in bioinformatics. It also covers important topics like RNA Seq Analysis, ChIP Seq Analysis in depth with code explanation of different real world examples which gives you head start in this beautiful subject.

The signature thing is frowned upon here too. I guess that is the reason why your post has been so heavily downvoted.

Thanks a lot.  Computational Genomics with R  is really interesting

One could read Bioinformatics: The machine learning approach as an overview of the established applications of machine learning to sequence alignment, RNA structure prediction, etc...
...Yet, this book is 20 years old, and there have been many important developments in bioinformatics and machine learning in this period, which have not yet been covered in any books. In other words, reading books on bioinformatics and/or machine learning is good for getting firm basic background in either of them, but any such books are inevitably out of date at the time of their publication.
I would like however to stress the importance of the basic background - one regularly sees announcements about breakthroughs achieved by reinventing techniques that were already described in the book that I mentioned.


text:
How to remove sequences from a fasta file using a sequence ID list which contains a space within the id?
I have a fasta file that contains sequence reads and sequence id file that needed to be removed from the fasta file. I have done this earlier, but since id contains a space my piece of code is not working and I don't know how to change it since I'm still learning. fasta file looks like this.
>m64041_200717_231916/74/12723_24868 id=30
CCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCA
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/105/20691_65844 id=488
AGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCAAGACTTTAT
>m64041_200717_231916/108/17414_67048 id=4956
TGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

IDList file is like this.
m64041_200717_231916/74/12723_24868 id=30
m64041_200717_231916/108/17414_67048 id=4956
m64041_200717_231916/105/20691_65844 id=488

I need to remove these sequences from that fasta file. I have a code I used earlier. But it is not working since id contains a space (I think). Previous code is like this.
# use as follows
# ./filter_seq.bash <fastafile> <list_of_ids_you_want_to_remove>
# ./filter_seq.bash all.fasta bacterial_IDlist

grep ">" $1 | sed 's/>//g' | sort | uniq > all_del
sort $2 | uniq > query_del
comm -23 all_del query_del | sort | uniq > filtered_del
perl -ne 'if(/^>(\S+)/){$c=$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' filtered_del $1 > filtered.fasta
rm *_del

I would be grateful if someone can help me.

I would highly recommend getting rid of the spaces first using something like sed 's, ,_,g' -I in.fa. Spaces in any kind of genetic data file are bad news and will cause downstream errors. It will make further analysis much easier.

@user438383 why do you say that? Spaces in fasta headers are incredibly common, allowed by the standard and should not cause any issues at all. I can't think of any tool that would have trouble parsing a fasta header like that and would even say that if a tool cannot do this, it is not fit for purpose. As an example, have a look at pretty much any fasta sequence from refseq. Here's human tp53: https://www.ncbi.nlm.nih.gov/nuccore/NG_017013.2?report=fasta.

Install (pip) biopython if u don't have it.
from Bio import SeqIO
with open('/path/ids.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip().split(' ')[0])
with open('/path/output.fa', 'w') as fout:
    for record in SeqIO.parse('/path/sequences.fa','fasta'):
        if record.id in to_keep:
            SeqIO.write(record, fout, 'fasta')

I tend to do this using the FastaToTbl and TblToFasta scripts I have posted before. FastaToTbl will convert fastq sequences to single-line, tab-separated recors which makes them far easier to parse. Using these two scripts, you can do:
$ FastaToTbl file.fa | grep -vwf idList | TblToFasta
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCC
TGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCAC
ACTTG

The -v option will tell grep to return non-matching lines, and the -f option tells it to read the patterns to search for from a file. In this case, we give it the file with the IDs to skip (idList). The -w ensures we only find "whole word" matches, so that something like m64041_200717_231916/77/1941_50622 id=31151111 will not count as a match for m64041_200717_231916/77/1941_50622 id=31151111.

Another approach would be to just use awk directly:
$ awk '(NR==FNR){s[">"$0]++} (/^>/){ s[$0] ? a=0 : a=1}a' idList file.fa
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

Here is the same awk script slightly expanded to make it easier to understand:
 awk '{
        if (NR==FNR){idsToSkip[">"$0]++}
        else if (/^>/){
            if($0 in idsToSkip){
                wantThisId = 0
            }
            else{
                wantThisId = 1
            }
        }
        if(wantThisId == 1){
            print
        }
      }' idList file.fa


text:
How to convert a SDF molecular structure file to a pandas dataframe with Python?
What is the best python tool to convert a SDF file to a structured data frame (pandas.DataFrame). Similar to Molconvert from ChemAxon? SDF can have different datafields for each compound which I would like to have in a dataframe. The 2D or 3D structures can be disregarded.

Please provide a small example file

It is easily done with RDKit:
from rdkit.Chem import PandasTools

fn = '/mypath/dataset.sdf'

df = PandasTools.LoadSDF(fn, embedProps=True, molColName=None, smilesName='smiles')


text:
blast p-error with making directory
I tried running this command python blastp.py CmydPro tlr TLR.txt to extract the query sequence for the target species I'm looking at. However, whenever I try running it always comes up as a blast p error and says there's no directory. I downloaded blastp.py as a file made from my lab's GitHub but I have no idea how to troubleshoot the error since it's not a required directory. I made sure all my directories were correct and that all my files were intact. Any help with this would be greatly appreciated.
This is the error I keep getting:
Traceback (most recent call last):
File "blastp.py", line 53, in
blastp_scout = subprocess.Popen(['blastp', '-db', organism_db, '-query', 'temp.fa', '-out', 'blast', '-num_descriptions', '5', '-num_alignments', '0'])
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 854, in __init__
self._execute_child(args, executable, preexec_fn, close_fds,
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 1702, in _execute_child
raise child_exception_type(errno_num, err_msg, err_filename)
FileNotFoundError: [Errno 2] No such file or directory: 'blastp'

These are all the files I have stored in my directory:
'$RECYCLE.BIN' 'Genome Files' TLR.txt
'ASM 28112' 'Genome Files (1)' Tham
'ASM 28964' 'Genome Files (2)' Vim
'ASM 3597' 'Genome Files (3)' Xenopus
'ASM 45574' 'Genome Files (4)' blast-2.11.0+
AnoCar Git blastp.py
Chrysemys GopAga get-pip.py
Cliv GopAga2.0-Genome human
CmydPro Halia ls
CmydPro_results.txt LICENCE mouse
Config.Msi Lib ncbi-blast-2.6.0+
CroPor Miniconda3-latest-Linux-x86_64.sh ncbi-blast-2.6.0+-win64.exe
'FASTAs (2)' OpenSSH-Win64 pip
Format Ornitho pvi
GCF_015237465.1_rCheMyd1.pri_genomic.fna PelSin pyzo
GFFs Python setup.py
'GFFs (2)' PythonX temp.fa
GRCh38 README.txt tlr_CmydPro.fa
GRCm38 Rnor tlr_ref.fasta
Gallus 'Rnor - Copy' vs_buildtools__576366542.1609221376.exe
GavGan Scripts vs_community__576366542.1609221376.exe
Gekko 'System Volume Information'

There is a windows version. It seems that I already downloaded blast but I need to do this now:
Command 'blastp' not found, but can be installed with: sudo apt install ncbi-blast+
Do you have any idea how to might download this?
I already downloaded blast as a tar file and the 2.11.0+ version. I'm just not sure how to download just ncbi-blast+ on its own.

That looks like you are using Windows. have you asked the developers if the script you are trying to run is even supposed to work on Windows? Is there a Windows version of blast? That can be launched from the command line?

Could you please tell us what operating system you're running (e.g. Windows, Mac, Linux)? If Linux, could you please show the output of the command uname -a?

You don't have blastp in your path. Blast is either not installed on your system or it's just not been added to path right.


text:
Trouble aligning next generation sequencing data to reference genomes using QuasR package in Bioconductor. Cannot import .txt
I'm trying to check the quality of my paired end read sequencing data. I am following this pipeline (https://f1000research.com/articles/4-1062#ref-21) which uses QuasR in the first step.
My list of paired files is stored in sampleFile. Currently this file cannot be read into the function. Not sure why. I'm using R studio on a windows operating system.
Also open to using other alignment or QC pathways if I can align to multiple different genomes.
auxFile <- "Ref_Genomes/auxiliaries.txt"
sampleFile <- "./HCV_rna_paired.txt"
genomeName <- "BSgenome.Hsapiens.UCSC.hg38"
td <- tempdir()

proj1<- qAlign(sampleFile, genome=genomeName, auxiliaryFile=auxFile)

Error: cannot open ./HCV_rna_paired.txt

Does the file named HCV_rna_paired.txt exist? Is it in the directory from which you run this R script?

Not a Windows user anymore for almost 15 years now but doesn't Win use backslashes to delimit paths and folders?

Yes the file exists in the same directory in which I was running commands.
As for the error message, would I receive the same message if it was a path error? I will try this next.

Definitely a file not found problem. Try: mydat<-scan("HCV_rna_paired.txt",what="",n=10) and see if it gives the same error


text:
How to retreive all the genes belonging to one GO term
I'd like to retreive all the genes belonging to a GO term.
I tried (in R):
library(org.Hs.eg.db)

xx <- as.list(org.Hs.egGO2ALLEGS)

xx is a list of GO terms, and for each GO term there's a list of the corresponding genes with an "entrezID" format.
The problem is that when I check this output on other databases (GeneCards, etc...) I can't find the same results and correspondence between the genes and the GO term that I found with org.Hs.eg.db.
Do you know a better way? Am I missing something here?

yes it's correct. how is it different from what you expect?

for example, the first GO term given by org.Hs.eg.db is "GO:0000002". The first gene in the list for this GO term is "142" which is "PARP1". PARP1 is not given as being part of this GO term in Amigo2 for instance or GeneCards.

you have to check what type of GO terms the websites are displaying, for example genecard will display only molecular function go terms and also certain evidences. This is not a programming problem. I suggest you post this on biostars if you still cannot resolve this

Which version of Bioconductor are you using? This will likely explain the difference. Also, not sure how frequent does GeneCard update their GO annotation or where do they take it from (It doesn't appear on the list of sources https://www.genecards.org/Guide/Sources)


text:
Is there a difference between any scRNA-seq visualization and pseudotime?
Is there a good definition of pseudotime? Some tools are clearly labeled as pseudotime and produce values along a trajectory, but there are more complex approaches that involve branching and disconnected populations. There are 70 of them in the dynverse, so there are many options.
UMAP will also order cells, but that is generally not considered pseudotime. You can even overlay Monocle or PAGA results onto a UMAP so the cells are still ordered the same way. SPRING is specifically designed to capture continuous trajectories, but is also not considered pseudotime. Why not?

Pseudotime is the order of cells along a trajectory which in the simplest case is a minimum-spanning tree based on a (low-dimensional) manifold such as PCA or UMAP. Without a defined trajectory you have no pseudotime, and therefore UMAP is not pseudotime, but in fact can be used to define the trajectory.

To simplify: a PCA/UMAP plot is not pseudotime, but if you draw a line over it, that would be pseudotime, right? All the cells are still in the same positions, so they still show the same relationships. Adding that extra coordinate makes it pseudotime?

Data scientists would argue that the "line" would need to come from a very elaborate machine learning approach, but yes, basically a line through the data makes it pseudotime. An accurate prediction of the root is critical though to give it some biological meaning, so the cluster of cells that most likely are the origin of a trajectory, e.g. stem cells.


text:
Retrieving all genes for a gene ontology term
I'm looking for an easy way to retrieve all the genes in a list that are associated with a certain GO term, preferably using R/Bioconductor packages. I'm not interested in under/overrepresentation or enrichment.
For instance, I want a list of all genes known to be located in 'presynaptic endosome' (GO:009883).
I tried the method referred to in an older post (Retrieve all genes for a specific gene ontology and https://www.biostars.org/p/52101/); following is my code:
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") #retrieve human.ensembl data via biomaRt

go_1 = c("GO:0098830", "GO:0098954", "GO:0099007", "GO:0099067", "GO:0099037", "GO:0098955", "GO:0099592", "GO:0099532") #GO IDs of all terms associated with "presynaptic endosome" (inclusive of child terms)

pre.gene.data <- getBM(attributes=c('hgnc_symbol', 'ensembl_gene_id', 'go_id', 'go_linkage_type'),
                    filters = 'go', values = list(go_1), mart = ensembl)

but my code doesn't filter the genes based on the GO terms. It gives me the following output:
GO:0006886 AP3D1 ENSG00000065000 intracellular protein transport IEA
GO:0016192 AP3D1 ENSG00000065000 vesicle-mediated transport IEA
GO:0030117 AP3D1 ENSG00000065000 membrane coat IEA

... and 49 other entries for the gene "AP3D1"

I'm unable to understand what might be the problem. I replaced the filter "go_id" to "go" because of the updated filters used:listFilters(ensembl).
Please help!

It is filtering by your GO terms, it is only giving you genes associated with those GO terms. However the GO terms in the output are not linked to the GO terms in your input, so it is giving you all the GO terms linked to the genes which are associated with the GO terms you filtered by.
Essentially BioMart is not very clever. It has a table of all the GO terms linked to the genes, it filters the genes to only include those linked to the GO terms you've given it. Then it takes that list of genes and queries that same table completely separately, getting the GO terms for each gene.

Thank you for the comment; It does make sense.


text:
Are there open-sourced graph-based variant callers?
I would like to compare GATK with graph-based variant callers. I have seen the Fast and accurate genomic analyses using genome graphs paper from SevenBridges and the paragraph tool by Illumina, though it does not look like folks use their platform or it is not open-sourced (in the case of SevenBridges's GRAF).
The idea is to see if graph-based approaches detect better both SNPs and structural variants compared with existing methods. Additionally, I believe that if there exists some novel graph-based aligner/assembler one could use the metadata downstream to identify variants.

Not sure about variant callers, but since you asked about aligners, you should check out vg.

@user438383 thanks, it's worth mentioning WhatsHap, GraphAligner, and Minigraph as well.

Yes, there are. There are some suggestions in the comments (VG, WhatsHap, GraphAligner, Minigraph).
However, to be clear, the current default variant calling algorithm used in the GATK (the haplotype caller) is a graph-based algorithm. It constructs a De Bruijn graph from mapped reads for local assembly around variants. Other tools that also follow this strategy include Platypus and freebayes.
You might also consider multi-color genome assemblers. For example, Cortex and Corticall use colored De Bruijn graphs to aggregate reads from multiple genomes and call variants.

I think there is a difference between using de Bruin graph and represent a genome as a graph. The colored graph assemblers are interesting. Thank you.

I completely agree.


text:
How to read E-value annotation on NCBI BLAST?
I am seeking clarification on how to read the  E value annotation on NCBI's BLAST results.
What I understand so far:

E -values are the expected value of the number of random querys of the same length as the input that would result in a even or better S-score on the same sized database. Since E-values are estimates, not probabilities , they can be lower than 0.

However, if NCBI BLAST retuns an e-value of say: "4e -19" do they mean:

$ {4 \cdot 10^{-19}} $ (this would not result in negative values )
or do they mean 4 E-19 as in : "four queries score better in S with a -19 difference"

It's a number in regular scientific format, so it's the former. 4e-19 = 4 * 10^-19

thanks for the reassurance! .. man, sometimes i stumble upon one source that uses a weird notation and i start to get quite confused.. The source of my confusion was Victor Flores comment on  https://www.researchgate.net/post/What_is_E-value_exactly_means_and_what_is_it_1e-63_representing

After rereading i could have seen earlier that he is apparently wrong.

He never says anything that resembles your second interpretation. How did you end up with that interpretation?

You are totally right, i confused "beyond" with "below" (not a native english speaker) . This should never have resulted in such a stupid confusion,  unfortunately the feedback at university is a little bit sparse due to the pandemic circumstances , so you had to bare with my undergrad incertitude. I honestly thank you for clearing it out nonetheless. On the upside i learned my lesson, took a step back and studied Altschuls formulation.

Don't worry about it, glad I could help. A lot of the long-winded description of these terms can make one confused but as you gain some practical experience, you'd understand these better.


text:
Supplementary aligments in VAF
This question has also been asked on Biostars
I have a doubt, are supplementary alignment usually considered when the variant allele frequency is calculated? Thanks a lot.
In my case I have some Illumina short data edited with Crisp-Cas9 and I want to calculate for a specific position all the frequency of the 4 nucleotides. I would like to understand if the supplementary reads aligned in that position must be taken into consideration.
I've used pysamstats to calculate and it consider them, but my doubt is more if it right to consider them, if it is normally done or not. We are doing it on human and we expect only a change of a single nucleotide variant in the region if the editing works.
I'm interested in what is done generally, if there is a "rule" to consider or not the supplementary alignment when the AF is calculated.
In my case I have some Illumina short data edited with Crisp-Cas9 and I want to calculate for a specific position all the frequency of the 4 nucleotides. I was trying to understand if the supplementary reads aligned in that position must be taken into consideration.
GATK HaplotypeCaller and Mutect2 take them into consideration? I'm not looking for a specific variant caller, but if I want to take one variant caller as "reference", I would probably look at them.
Thanks a lot.

It depends on the tool and the use case but generally they should be.
It's not as important with short reads as only a small portion of reads have supplementary mapping (unless you're looking at breast cancer post chromothripsis, for example). But imagine the case where you're working with long reads (as I do). I have samples with a mean read length of 60kb. The human reference is structurally different to any given sample meaning that the chance of any given long read having supplementary mapping is high. If supplementary mapping didn't count towards allele counts then many would be artificially low. This is my use case, others are different.
I would count them unless you have some reason to suspect that they are are artefact of the editing and or the mapping. Sorry to never give a straight answer but that's bioinformatics! From what you described I would guess 99% are real but there's normally some noise and if that's of concern your code should include some QC. I.e, consider if multiple reads support the sup mapping and their quality.

Thanks you so much Liam!


text:
How to perform bowtie2 analysis with slurm?
I am trying to run my alignment script that works locally, using SBATCH.
Official manual for bowtie2 says I can use samtools for that. However when I run the following script I have an error
line 22: bowtie2: command not found

Below is the script:
#!/bin/bash

#SBATCH --time=05:00:00
#SBATCH --account=def-myaccount
#SBATCH --mem=200000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load samtools/1.9

for FILE in *R1_001.fastq.gz;

do
        BASENAME=${file%%_R1_001.fastq.gz*}

        NUM=`echo $FILE | sed 's/_S.*$//'`

        bowtie2 -x cov_nastya -U $FILE -S ${NUM}exit.sam

done

It seems you do not have bowtie installed or loaded. Use module load bowtie2 together when you load samtools.

Thank you, it worked. Just a general inquiry at the end.. Is samtools really necessary here? Thanks.

From your script samtools does not seem to be necessary.


text:
How can I use my Myheritage DNA results file for further analysis?
I had my DNA tested by Myheritage and they sent me a csv file with RSID, Chromosome, position and result (which base) with about 700,000 rows.
I understand most analyses of DNA use VCF files but is there anything i can do with this csv file i.e. check for genetic health-related genes?
I am not a bioinformaticist but I am a scientist and I can use R and python.  Ive heard of the Gnomad database but not sure if i can match things to my csv file.

I think people often use promethease for this kind of thing.

Check out data.table and merge with gnomAD or better, CLINVAR. That might give you some pointers.

A csv file with RSID, chromosome, position and result is enough for what you want to do and these are the core columns of a VCF (which is just a TSV with some headers describing how it was made).
Given you have ~700k rows I suspect that yes there will be genetic health-related SNVs (single nucleotide variant). Disclaimer; these data could include info pertaining to your health and your family's health, I strongly advise that you speak to a genetic counsellor to understand these things.
Gnomad is a population database so generally speaking (plenty of exceptions of course) if there is a SNV in there more than a few times its probably not causing disease. https://www.ncbi.nlm.nih.gov/clinvar/ and https://omim.org/ are example of disease databases.
Here's a pandas example with clinvar (https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/ you'll need to match your version).
import pandas as pd
df = pd.read_csv(
    'clinvar.vcf.gz',
    sep='\t',
    comment = '#',
    header=None,
    names = ['CHROM','POS','ID','REF','ALT','QUAL','FILTER','INFO']
)
df.head()

CHROM   POS ID  REF ALT QUAL    FILTER  INFO
0   1   930188  846933  G   A   .   .   ALLELEID=824438;CLNDISDB=MedGen:CN517202;CLNDN...
1   1   930203  972363  C   T   .   .   ALLELEID=959431;CLNDISDB=MedGen:CN517202;CLNDN...
2   1   930248  789256  G   A   .   .   AF_ESP=0.00347;AF_EXAC=0.00622;AF_TGP=0.00280;...
3   1   930275  969662  T   G   .   .   ALLELEID=959432;CLNDISDB=MedGen:CN517202;CLNDN...
4   1   930336  843786  G   A   .   .   ALLELEID=824439;CLNDISDB=MedGen:CN517202;CLNDN...

Then lets say one of your 700k SNVs is at chrom 1, position 930336 and you have an A
print(*df.query('(CHROM == 1) & (POS == 930336) & (ALT == "A")')['INFO'].str.split(';'))

Gives a list from which you could pull out CLNSIG which here is Uncertain_significance
['ALLELEID=824439', 'CLNDISDB=MedGen:CN517202', 'CLNDN=not_provided', 'CLNHGVS=NC_000001.11:g.930336G>A', 'CLNREVSTAT=criteria_provided,_single_submitter', 'CLNSIG=Uncertain_significance', 'CLNVC=single_nucleotide_variant', 'CLNVCSO=SO:0001483', 'GENEINFO=SAMD11:148398', 'MC=SO:0001583|missense_variant', 'ORIGIN=1']

So you could easily parse your 700k SNVs and get their respective clinical info this way. Please think about this carefully before you do it though, there are many implications including life insurance and mental health! I'm of the opinion that this sort of data is useful but MANY people disagree. Further, this advice is general in nature only, I accept no liabilty for what you or anyone else chooses to do with this publicly available data.


text:
PDB: Download all biological assemblies
I was wondering if it is possible to download all biological assemblies from the PDB at once (so not the standard asymmetrical files). For example, I want a file containing the coordinates of all atoms in a trimer, instead of only the coordinates of a monomer.
I saw that you can batch download with a shell script, but I'm not sure if you can also do that for biological assemblies.

Alternatively, you could download a "standard" file from the PDB (either in the PDB format or in mmCIF) and generate assemblies locally with one of many programs/toolkits have such capability.

I'm not aware of a program/toolkit with which you can do that automatically for a large set of PDB files. If you know of a way to do that then that would be perfect.

To process multiple files you will probably need a simple script. Note that PDB entries often define several alternative biological assemblies. As for toolkits, I'm most familiar with gemmi (because I maintain it): you could generate assemblies from the command-line with gemmi convert --assembly=1 ... or from Python. I suppose that BioJava, cctbx, PyMOL, etc also have such a capability.

Thanks a lot, I'll give that a try.

If you go to the wwPDB ftp site https://ftp.wwpdb.org/ in your browser you can find the biological unit structures (and anything else available in the PDB). The caveat of course is that a biounit structure has to have been generated for your structure.
By way of some examples, to get all biological units you could do (note the ftp instead of https in front of the URL):
wget ftp://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/all/*

or to get a version organized in directories:
wget ftp://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/*

Note that this would take a few days to download everything, though. After the first download, you can add a -N to the wget command to just download new files.
For getting just one file:
wget ftp://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/00/100d.pdb1.gz

Hope this helps.


text:
Obtaining Whole Genetic Sequence
As an end-user of my own data, specifically raw DNA sequence (WGS, whole genome sequencing). How or where do one obtain such DNA data so that a biology hobbyist can perform bioinformatics on this?
What keywords is best describe to the search engines when finding these online WGS websites).
It would be nice to have a USB stick mailed to me after such DNA HomeKit.  What file format is most commonly offered? (VCF?)
Am I looking too far ahead into the future?

Not sure if I follow, are you asking how or where you can get sequenced and what data you can expect when you would order such service?

something that a home biology hobbyist can work with.

It is not clear to me where you wanna start. Do you have a biological sample or already extracted the DNA? What is your source organism? Is it a prokaryote or eukaroyte? What kind of bioinformatic task do you plan to do?

Are you looking too far ahead into the future? No. This is certainly possible now. The gold standard being Illumina or BGI short read WGS. Long read sequencing can capture some extra data but is very noisy.
How or where do one obtain such DNA data? There are now many companies. You want one the uses NGS with reads at least 150 bp and garantees at least a mean of 30x read depth. In my opinion you also want one that gives you ownership of your data. Ideally they will give you; fastqs (raw reads), a BAM (mapped fastqs) and a VCF (variant file). I good company will ask you if you want to donate your data to science. Be warned these are not defaults. Many (most?) companies won't give you your raw data citing medico legal risk and then they will on sell it to a pharma companies. I think these guys are OK (as they are in Europe with good privacy laws) https://www.dantelabs.com/ but plz do your own research.
Last, these data could include info pertaining to your health and your family's health. Always good to speak to a genetic counsellor to understand these things.

My ownership is a big plus

"DNA data" comes in several forms dependent on the technology used to produce them.
Companies like 23andMe are using SNP chips and those are available for an only rather a limited number of species. To be completely honest, I don't even know what would be the raw data they would provide, but I would suspect it could be a vcf file.
If the chosen technology would be DNA sequencing. DNA sequencing is provided by a few dozen of sequencing centres around the world, some of them associated with universities. DNA sequencing generates files that contain individual sequencing reads with assigned quality scores ("how much the sequencer trusts each base"), the typical formats are fastq (Illumina) or .bam (long read technologies). Usually, you would also need to find a suitable reference in one of the public databases, for instance, the NCBI portal for human reference. Then map your reads and do variant calling. Of course, then it depends on what would you be interested in.
If you don't care for your own DNA and you are more interested in life around you, you can check ncbi for all sorts of genomic data. If you would be interested to understand how the raw sequencing data work, ebi is a good place to look. Also, in every single paper generating/using sequencing data you there should be mention where in public databases you can find them.
Finally, let me warn you. Some of the genomic analyses are quite computationally heavy (not all though). If you don't have a specific question in mind, it would be wise to start with bacterial sequencing.
Hope you will have fun with DNA data!


text:
How can I run the latest version of samtools on Slurm?
I am trying to run the following script on slurm to extract information from .sam files. While everything is perfect on my local machine with the samtools version 1.11, it does not work on slurm. I have this error:
samtools: command not found

That means slurm cannot understand the module if I am not mistaken. What would be a solution here?
#!/bin/bash

#SBATCH --time=05:00:00
#SBATCH --account=def-myaccount
#SBATCH --mem=200000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
# module load samtools/1.11 # this is what I tried first
module spider samtools/1.11 # this does not help either...

for FILE in *.sam;

do

        NUM=`echo $FILE | sed -r 's/.sam//' | sed -r 's/nastya//'`

#         echo "$FILE" ;

        samtools flagstat -O tsv $FILE > file${NUM}.txt

done

@OP, these kinds of specific questions should be discussed with your sysadmin as it is specific to the cluster you work on. Alternatively you can install a local version with package managers such as conda but.

I would start with:
module avail

Then look to see (a) if there is a samtools module avilable and (b) exactly what it's called.
If it's called samtools/1.11 then:
module load samtools/1.11

Should work, if it doesn't speak to your sys admin.

Great. Thanks for the tips. It worked. Didn't know about avail.


text:
Multiple box plot
I am using below code to plot multiple box plot but I am still struggling to understand how I can add whisker to box plot and how to change sample order like in this figure SampleType?
library(reshape2)
library(ggplot2)
list.files()
morphData <- read.table("Box.plot.chemical.txt", sep="\t", header=TRUE)
morphData
df.m <- melt(morphData, id.var = "SampleType")

pdf("box.plot.pdf", width = 18, height = 12)

p <- ggplot(data = df.m, aes(x="", y=value)) +

theme_bw() +
theme(
      panel.grid.major=element_blank(),
      panel.grid.minor=element_blank(),
      plot.title = element_text(vjust = -8.5,hjust = 0.1),
      line=element_line(size=1),

      axis.text.x = element_text(size = 10, face ="bold", hjust = 0.5, colour = "black"),
      axis.text.y = element_text(size = 12, face ="bold", hjust = 0.5, colour = "black"),
      axis.ticks.x = element_blank()) +

theme(strip.background=element_rect(fill="white")) +
theme(strip.text=element_text(color="black",face="bold", size =14)) +
geom_boxplot(aes(fill=SampleType),outlier.colour = NA)
p + facet_wrap(.~variable, scales="free")
dev.off();

update:
I was able to remove the variable name from the bottom by using above updated code but still getting tick mark on x-axis, and also how I can make these label bold?
update:
It will be a great help if you could also guide me how I can add Tukey's posthoc significance letter in above box plot to make this plot more meaningful.  I was able to make bold facet wrap text, now please help me with posthoc test letter. I have updated the code again.
update:
Can you please help to change the order of samples; like I want to change the order of SampleType in this box plot and also how we can remove outlier and will it be possible to add whisker so it will look better. I was able to remove outlier by using outlier.colour = NA, Now It just whisker, but I am not able to understand How I can add whisker to box plot? I have updated the code and figure above.
update:
sorry, I wanted to add errorbar and with this command, it is not adding error bar on all box: stat_boxplot(geom ="errorbar")

This is a good question and the answer below is very useful. Please could you consider accepting the answer?

what do you mean by How I can put the variable name on the top of each box plot like Rice_P at top and at the bottom as well? Do you want to keep the label at both top and bottom or remove one of them? From my opinion it's not necessary to use facet in this case, since you only have one variable on x-axis in each facet now. You can just put variable as x-axis, and color your box as it is now. Boxes with the same variable are going to be grouped together.
You can use axis.ticks.x in theme to remove the ticks on x-axis similar to axis.text. To make label bold maybe you can try fontface="bold"
Adding testing letter may be harder, you could check out ggpubr or ggprism package.
you can use as.factor() to change sample name to factors and plot them in order. Or you can use scale_x_discrete(limit = <your sample order>) to force an order. I don't see outlier in plot. Outliers are usually black dots in box plot. And the box plot you show already has whiskers, so I 'm not sure what you want to add.
I'm not sure whether it's a good idea to plot error bars in boxplot because that is not the convention. But probably you can try using geom_errorbar in addition to geom_boxplot as another layer but it requires more customization.


text:
How can I use pandas agg here to avoid iteration?
can anyone help me improve this pandas code?
import pandas as pd

df = pd.DataFrame(
    [
        [
            'chr1', 222
        ],
        [
            'chr1', 233
        ],
        [
            'chr1', 2123
        ],
        [
            'chr2', 244
        ]
    ], columns = ['chrom', 'pos']
)

df2 = pd.DataFrame(
    [
        [
            'chr1', 221, 223
        ],
        [
            'chr1', 230, 240
        ],
    ], columns = ['chrom', 'start', 'end']
)

Gives me 2 dfs with genomic coordinates. The first one is an exact position:
    chrom   pos
0   chr1    222
1   chr1    233
2   chr1    2123
3   chr2    244

and the second is a range:
    chrom   start   end
0   chr1    221 223
1   chr1    230 240

I need to find the count of exact coordinates that are in one of the ranges (in the same chrom)
This works but is slow:
c=0
for chrom, data in df.groupby('chrom'):
    tmp = df2.query(f'chrom == "{chrom}"')
    for p in data.pos:
        for s, e in zip(tmp.start, tmp.end):
            if s < p < e:
                c+=1

Then c = 2
I think I can use agg to do this without iteration (and hopefully faster) but I can't get it working. Can anyone show me how?
PS I am also asking this on stackoverflow.

https://codereview.stackexchange.com/ is probably the best place to ask questions about optimising code.

That actually a really good question. A bioinformatics question. It's self contained and we'll explained. If U don't know just don't answer. No need for negatively.

I'm actually trying to help you, since code review (hint) is a much better place to get your code reviewed than either here or SO. You seem to be new here - so different sites have different functions. It's not about negativity / whether or not I can answer the question, it's about the correct questions are on the correct sites.

Ok fair enough.

# pip install pyranges or conda install -c bioconda pyranges

import pyranges as pr

g1 = pr.from_string("""Chromosome Start End
chr1   222 223
chr1   223 224
chr1   233 234
chr1   235 236
chr1   2237 238
chr1  2123 2124
chr2   244 245""")

g2 = pr.from_string("""Chromosome Start End
chr1    221  223
chr1    230  240
chr2    0 1000""")

r = g2.count_overlaps(g1)
r
# +--------------+-----------+-----------+------------------+
# | Chromosome   |     Start |       End |   NumberOverlaps |
# | (category)   |   (int32) |   (int32) |          (int64) |
# |--------------+-----------+-----------+------------------|
# | chr1         |       221 |       223 |                1 |
# | chr1         |       230 |       240 |                2 |
# | chr2         |         0 |      1000 |                1 |
# +--------------+-----------+-----------+------------------+
# Unstranded PyRanges object has 3 rows and 4 columns from 2 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

r.df
# 0       chr1    221   223               1
# 1       chr1    230   240               2
# 2       chr2      0  1000               1

How did I not know about pyranges!? Thanks!


text:
co-occuring genes in my list
I have a files which has columns representing genes,sample numbers and sample names (55). The genes are displayed by occurrence in the samples (highest to lowest-55 and 2).I wish to calculate the genes which may co-occur in samples.
I tried doing this with awk 'n=x[$1,$2,$3,$4,$5,$6]{print n"\n"$0;} {x[$1,$2,$3,$4,$5,$6]=$0;}' adding from 1 to 55 columns, but it doesn't work.
An example here:
Nfib    55  A1.S1   A2.S2   A3.S3   A4.S4   A6.S6   A7.S7   A8.S8   A9.S9   A9.S9*  A10.S10 A11.S11 A12.S12 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.8   IZ_BH_2.9   IZ_BH_2.10  IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.1   IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.7   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_3.11* IZ_BH_3.12  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.4   IZ_BH_4.5   IZ_BH_4.6   IZ_BH_4.7   IZ_BH_4.8   IZ_BH_4.9   IZ_BH_4.10  IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.2   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

Lyn 41  A1.S1   A4.S4   A7.S7   A8.S8   A9.S9   A9.S9*  A10.S10 A11.S11 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.9   IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.3   IZ_BH_4.4   IZ_BH_4.7   IZ_BH_4.8   IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

Cdk5rap2    52  A1.S1   A2.S2   A3.S3   A6.S6   A7.S7   A8.S8   A9.S9   A9.S9*  A11.S11 A12.S12 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.8   IZ_BH_2.9   IZ_BH_2.10  IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.1   IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.7   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_3.11* IZ_BH_3.12  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.3   IZ_BH_4.4   IZ_BH_4.6   IZ_BH_4.8   IZ_BH_4.9   IZ_BH_4.10  IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.2   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

The first name is of the gene, followed by number of samples and sample names.
for these 3 genes these are the commonly occuring samples:
A1.S1, A7.S7, A8.S8, A9.S9*, A9.S9, A11.S11, IZ_BH_2.1, IZ_BH_2.2, IZ_BH_2.3,IZ_BH_2.4, IZ_BH_2.5,IZ_BH_2.6,IZ_BH_2.7, IZ_BH_2.9, IZ_BH_2.11,IZ_BH_2.12, IZ_BH_3.2,IZ_BH_3.3,IZ_BH_3.4, IZ_BH_3.5, IZ_BH_3.6,IZ_BH_3.8,
IZ_BH_3.9,IZ_BH_3.10,IZ_BH_3.11,IZ_BH_4.1,IZ_BH_4.2,IZ_BH_4.4,IZ_BH_4.8,IZ_BH_4.11,IZ_BH_5.1,IZ_BH_5.3,IZ_BH_5.4,IZ_BH_5.5,IZ_BH_5.6,IZ_BH_5.7,IZ_BH_5.8

Kindly help.

Show a small worked example plz

I have added an example

I think what makes this especially hard is the way your data is formatted.
The problem would be way easier if you had a matrix of presence / absence instead. Each row would be a gene, each column a sample, and the values would be 0 (absent) or 1 (present).
Awk becomes a bit inconvenient for this task, so I'll use python instead. To get such a matrix from your file, you could do the following:
import numpy as np
import pandas as pd

# Get total number of genes / samples
all_samples = set()
all_genes = set()
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    all_genes.add(line[0])
    all_samples.update(line[2:])
all_samples = list(all_samples)
all_genes = list(all_genes)

# Initialize presence matrix with correct dimensions
pres_mat = np.zeros((len(all_genes), len(all_samples)), dtype=int)

# Map each sample/gene to a row/col
row_map = {gene: i for i, gene in enumerate(all_genes)}
col_map = {sample: i for i, sample in enumerate(all_samples)}

# Fill matrix incrementally
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    curr_gene = line[0]
    row = row_map[curr_gene]
    for sample in line[2:]:
      col = col_map[sample]
      pres_mat[row, col] = 1

# Add names and write to file
df = pd.DataFrame(pres_mat, columns=all_samples)
df.insert(0, 'gene', all_genes)
df.to_csv('gene_presence.tsv, sep='\t')

You should get something that looks like:
       gene  IZ_BH_5.5  A8.S8  IZ_BH_2.1  A9.S9*  IZ_BH_5.3  A9.S9  ...  IZ_BH_3.11  IZ_BH_4.7  IZ_BH_4.9  IZ_BH_2.8  IZ_BH_3.5  IZ_BH_2.6  IZ_BH_4.3
0       Lyn          1      1          1       1          1      1  ...           1          1          0          0          1          1          1
1      Nfib          1      1          1       1          1      1  ...           1          1          1          1          1          1          0
2  Cdk5rap2          1      1          1       1          1      1  ...           1          0          1          1          1          1          1

With this format, it is now easy to compute the most frequent samples. E.g. df.sum(axis=0) to get the number of occurence of each sample for these 3 genes. It would also allow you to do different things, such as computing the pairwise similarity between samples.

Thanks a lot @cmdoret


text:
Do RepeatModeler results contain functional domains?
The repeat families predicted by RepeatModeler contain known transposable elements (TEs) and unknown ones.
How do we know whether some of these may actually:

within a functional domain of a gene or
represent a functional domain?

Could it be predicted because repeats may simply due to this domain appeared in multiple gene loci?
Thanks.


text:
How to determine a correct promoter sequence? (Bisulfite sequencing)
I'd like to ask you for help with an university assignment. First let me post the entire question:

Design primers for bisulfite sequencing of the promoter of the ACE2
gene in microbat  Myotis lucifugus (based on the annotation of the
Broad Institute Myoluc2.0/myoLuc2). The gene is localized on the
contig GL429816. As a transcription start, use the start of the ACE2
transcript variant XM_023753670.1; once you find ACE2 in the UCSC
genome browser .You will find it in NCBI RefSeq track

I don't want anyone to solve it for me but I do have a problem with one aspect. I am not clear on how to determine what is a valid promoter for this gene in microbat. I do know it is the upstream sequence however I am not clear on how to tell the number of upstream bases that are the actual promoter. I did find some sources for the actual promoter sequence determined by much smarter people than myself however those are for ACE in humans so I am afraid that it is not applicable in my case.
I'd like someone to point me in the right direction where I could find my answer would be enough.

Thank you I decided to use 900bp, just because that is the legnth the teacher used in of his lectures. I would like to select you comment as an answer, but it's not possible it seems.

Posted it as answer, not sure if a moving function exists, maybe if you are a moderator.

I think there is no fixed size. Ideally you would have something like histone ChIP-seq or ATAC-seq data to approximate the relevant upstream interval. Since this is not always possible most people go with -500bp. Others might have different suggestions though. Maybe matching known motifs of transcription factors that regulate that promoter to the upstream region can help approximating a better window. But as Sanger sequencing usually gives about 1kb of sequence I'd just take a 500bp to 1kb window and sequence that, fwd and rev primers for sequencing to have double confirmation.


text:
Not Able to pull genes due to file.fa error
I've recently been pulling TLR genes for different species and keep running across an error called file.fa whenever I start pulling them. Whenever I pull the first gene it's fine and makes the file.fa but when I pull the second gene (TLR4) it comes up with a file.fa error and deleted the previous file.fa
made from when I pulled the first gene. Any feedback would be greatly appreciated.
['NP_001273862', 'XP_006138208', 'XP_006138164', 'XP_014427015', 'XP_014429739']
Traceback (most recent call last):
  File "blastp.py", line 104, in <module>
    SeqIO.write(record, "file.fa", 'fasta')
  File "/home/dylan_w190/miniconda3/lib/python3.8/site-packages/Bio/SeqIO/__init__.py", line 522, in write
    with as_handle(handle, "w") as fp:
  File "/home/dylan_w190/miniconda3/lib/python3.8/contextlib.py", line 113, in __enter__
    return next(self.gen)
  File "/home/dylan_w190/miniconda3/lib/python3.8/site-packages/Bio/File.py", line 72, in as_handle
    with open(handleish, mode, **kwargs) as fp:
FileNotFoundError: [Errno 2] No such file or directory: 'file.fa'

Can U plz provide enough detail for me to try to reproduce the error?

Hi there, you seem to be new to this forum. It is good practice to provide full code to actually reproduce this error. Otherwise your question will quickly get closed. Right now we neither know about your input files, not the code you used and whether this is a tool or custom script. Please add all that.
This here should be the starting point to investigate:
FileNotFoundError: [Errno 2] No such file or directory: 'file.fa'

Show your code please- the error is slightly counterintuitive to me since writing to a file should not lead to FileNotFoundError


text:
Do you know any virtual screening libraries with small, soluble compounds?
I had a quick question on virtual libraries. I'm hoping to perform a virtual screen using a library that contains very small, soluble compounds - think glycerol and smaller (<100MW). As far as I can tell, the majority of virtual libraries contain mostly larger, drug-like compounds. Is anyone aware of any libraries that might contain ideally 10,000+ or so of these tiny soluble compounds?
If not, any idea what the process might be (and how difficult it would be) to construct a custom virtual library enriched for these compounds?
Thank you!!

Hi, how is the htslib tab relevant here?

Lipinski and logP
Most datasets will be filtered by the Lipinski rule, plus in the some cases, such as Enamine Real whether they can be bought.
One criterion in the Lipinski rule filter is the logP, or octanol-water partition coefficient. This is super important, because a compound that cannot cross the membrane is likely of limited uses —as a consequence prodrugs are a commonly taken approach.
However, albeit obvious, not all datasets in Zinc are filtered by the Lipinski rule, they just need more browsing...
Natural?
However, you may have a valid reason to eschew logP.
Nature tries its best to have intermediate metabolites that have hydrophilic handles attaches, eg. CoA, phosphates and sugars. Therefore, if you want a natural metabolite this is a completely different question. Many of these are not very purchasable, but a PDB-derived dataset of natural ligands will likely be purchasable or easily extractable. Alternatively Zinc has several natural catalogues...
Too small
Glycerol is a 3 carbon molecule and has a MW of 90 Da. Many of the old drugs are small, but not that small. For example aspirin is 180 Da, which is twice the size of glycerol. Lipinski rule says not to go over 500 Da and 250-400 Da is a lovely range. Say you have a perfectly bound glycerol, with a perfect ligand efficiency of -1.5 kcal/mol/atom, that would mean it would bind at –9 kcal/mol in the best case scenario, i.e. it is too small to bind anything with nanomolar affinity and will actually bind a lot of things. Glycerol after all is commonly cocrystallised in protein...
The GDB datasets have enumerated compounds of a given size, but don't go that low.
Docking
There will be two further complications that will hinder the accuracy.

protonated and deprotonated atoms
polarisable groups

The former requires correct conformers, the latter requires a force field that has a Drude particle, e.g. in openMM.
Caveat against sorting by binding score
Also, it should be said that if you sort by binding score (kcal/mol or kJ/mol) you will get bigger compounds (hence the whole point of ligand efficiency) and compounds with high logP. The former is because they will have more contact points to make, the latter is because given that ∆∆G is ∆G_bound - ∆G_unbound, an unbound high-logP compound is highly unhappy in aqueous solution.

ZINC20 is a free database of commercially-available compounds for virtual screening. It contains 1.3M compounds with molecular weight less than 200.
https://zinc20.docking.org/
Cartblanche22 is an interface to ZINC22. Here you can find more compounds than ZINC20.
https://cartblanche22.docking.org/cartblanche


text:
multi--platform GEO analysis based in expression of single gene
This question has also been asked on Biostars
I have a challenge in analyzing GSE34941( It has 4 platforms) according to the expression of a specific gene in R. First, which is the best way to group my patients according to the expression of a gene to high- and low-expression? Second, which is the best way to remove the batch effect of multi-platform analysis here?
I really appreciate any help you can provide.

Maybe can try the ComBat to remove the batch effect


text:
Differential expression analysis for a subset of transcripts
I want to check the differential expression of a specific class of transcripts (say, long non-coding RNAs) using DESeq2. Now, I know that the normalization step takes into account the total number of reads per sample. If I only input a subset of all the transcripts, the sequencing depth will change.
So my question is: can I simply input the subset of interest? Or should I input all transcripts, and subset the results later?

It's safer to import everything.  You'll want all the data for normalization and dispersion estimates.