text:
Find a genomic coordinates for a protein aminoacid position
Is there a function that can map genomic position (hg19) back from a protein position? I can have name of a particular transcript, and exon number.
For example, I have KRAS gene for which I would like to know the exact genomic position of G12, and I expect to have an output saying:
chr12, from 25398283 to 25398285.
transvar could be an excellent solution if the current built was working (broken links for fasta annotations). The web search is working but not allowing files heavier than 5kB which is my case.

This is considerably more complicated since there can be more than one genomic position referring to G12 in the protein in cases where multiple transcripts use different exons and yet >1 of them result in a protein with a G at position 12.

@lizaveta could you please review the answer? In my opinion this does answer your question (just my opinion), if it does please kindly consider "accepting" the answer.

You can do this manually by looking up the protein variant on https://varsome.com. For your example, that would mean searching for KRAS G12V (change V for anything else, the idea is to give an HGVS protein notation for a variant, so any single letter amino acid code would work). That would take you to this page https://varsome.com/variant/hg19/KRAS%20G12V?annotation-mode=germline, which includes:

Of course, I imagine you want to automate this. VarSome can help there too, through our API, but that is a paid service. If you have the option of paying for this, the VarSome API can be used to (among other things) look up arbitrary variants and translate between different notations, so you could easily adapt this to get the information you need.
Disclosure: I work for the company behind VarSome.


text:
Gene not found in Affymetrix expression profiles
I am studying the ABA network in A. thaliana, consisting of HB7, ABI1 and AREB2. The AGI code I was given are, respectively: AT2G46680, AT4G26080 and AT1G45249.
I downloaded the following file in order to convert the array element name to the AGI one, but I am not able to find AT1G45249 in the above file. How can I check whether AT1G45249 has another name? Or even better, how can I check what is the array element name in Affymetrix (GPL198) corresponding to AT1G45249. Thank you.

You can use  org.At.tair.db: Genome wide annotation for Arabidopsis, which is an R package available on Bioconductor to extract gene symbols.

source: Stack Exchange

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

sort: votes


text:
Converting Plink .bed files to .bcf *with* an AC (allele count) field?
I have some extremely large (N ~ 400,000) Plink .bed files that I want to phase using shapeit5.
Shapeit5 requires input data be in .bcf or .vcf format, and that it have an "AC" (allele count) field.
Here is what I have managed to figure out:

I can convert the files to .bcf easily enough with Plink2's --export bcf option, but the resulting .bcf file has no "AC" field.
I can add an "AC" field with bcftools +fill-AN-AC or bcftools +fill-tags -t AN,AC which outputs a new file with that field added, but this is extremely slow (several hours) because bcftools only uses one thread to calculate the allele counts (the --threads flag is for output compression only, it does not affect these commands)
I can split the chromosome into several chunks, use bcftools to add an "AC" field to each chunk, and then use bcftools concat to combine the chunks back into one big .bcf file. This is much faster, but also undesirable because of the amount of storage space used (requires temporarily making 3 copies of a very large dataset)

Is there a tool that will let me go from Plink .bed to .bcf while adding an "AC" field all in one step? Or if not, is there at least a better option for adding an "AC" field to the large .bcf file that will run with multiple threads and not require me to manually split the data into chunks?

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


text:
How drastic would the result be if I use fasta or reference assembly from ucsc and gtf from gencode?
There are difference annotation file for UCSC and gencode.
But if I use the reference assembly from UCSC and the GTF from Genocode or vice versa would my downstream results would be wrong?

Do you want to include older UCSC/gencode releases or only the most recent ones? They should now be fully compatible.

so i used gencode.v21 as my gtf file and not sure about the reference assembly,where i have used hg38 but im not sure its from gencode or ucsc .Is there a way to find out the what reference assembly im using from the reference assembly ?

Never use genomes or annotations from UCSC, they're poorly versioned and only recently with mouse and human have they even included all of the contigs. For fasta/GTF files from early in the GRCh38 release, you can tell whether you're using UCSC or Gencode by the presence/absence of _random contigs, which will only exist for UCSC. These were mostly later split into the actual contigs, so recent download from UCSC should more closely match what you find at Gencode/Ensembl. Further, that time predated UCSC beginning to adopt Gencode's vastly superior annotations, so if your GTF file has instances where the same gene ID is on either multiple strands or multiple chromosomes (this is obviously biologically impossible) then you have a UCSC GTF file.
In general, with early GRCh38 releases your only real issues will be with _random contigs, which are a minority of the genome and don't have all that many genes. But really, you should be keeping track of the sources of your files and ensuring that they're compatible.
Update: I should expand a bit on my "Never use genomes or annotations from UCSC" comment. In point of fact the genomes themselves aren't so terrible. Early on UCSC had the bad habit of concatenating contigs together into _random "chromosomes", but they seem to have mostly kicked that habit as of late. Note, however, that there are no versions for their genomes. Since reference genomes continue to get updates over time (mostly through the addition of patches) the lack of actual versions means you have to manually check if a recently downloaded file matches what may have been downloaded either previously or by someone else. This has obvious consequences when it comes to reproducibility. The same issue occurs for annotations from UCSC, but they have the additional problem of historically having biologically incoherent concepts of genes. That is, they will contain the same gene in multiple places with multiple orientations, which will break many tools in both obvious and completely unclear ways. For example, DEXSeq will simply break with an error message if given a UCSC annotation, since they break biological plausibility. If you were to use these annotations files with deepTools, you wouldn't get an error message, but the resulting output would be only partial, due to the biologically impossible annotation effectively corrupting most obvious ways of storing annotation data in a data structure (i.e., you can no longer treat IDs as unique). This could have downstream ramifications on biological interpretation of results.

to make sure if im have gencode so i used this "cat gencode.v21.annotation.gtf | grep "_random" no instance of that word .

You could also just head the file. Gencode files tend to start with a few comment lines (granted, given what the file is named, the odds of it being from Gencode were very high to begin with).

##description: evidence-based annotation of the human genome (GRCh38), version 21 (Ensembl 77)
##provider: GENCODE
##contact: gencode@sanger.ac.uk
##format: gtf
##date: 2014-09-29   im sure now im in right path

@krushnachChandra not really relevant, but just so you know, grep can read files. You don't need cat file | grep foo, you can always do grep foo file.

sure i was just trying to see if i haven;t messed up my annotation file .I will surely keep that in my mind. "grep"

Devon, could you elaborate on your advice never to use UCSC genomes. Never seems a bit extreme. Can you explain what they're missing in a bit more detail?Is there any reason not to use them for human reference? Yes, they may have *random in them, but is that such a problem?

*random is the least of the problems. I'll elaborate more in my answer.

Thanks, that makes sense. I have used the (human) UCSC genome sequences extensively and haven't had any issue with them although I admit that the lack of versioning is a real pain. I was afraid you were suggesting the actual sequences were somehow of lesser quality.

The main difference is in the way the chromosomes are named - UCSC uses the "chr" prefix (so chromosome 1 is "chr1") while in gencode the "chr" isn't used (so chromosome 1 is just "1"). Depending on your use case, this can obviously cause problems - if you're trying to match a locus (e.g. from gencode 1:1000002) between them, whatever tool you use is going to be looking in your aligned data for "1:1000002", but in your aligned data, it'll be named "chr1:1000002", so it won't match the two up.

Both use the chr prefix for the most recent human and mouse releases.

so its with the annotation file my results would vary  but not much with the references?


text:
TCGA data sets joining multiple files
I have downloaded harmonised the TCGA-HNSC dataset and extracted it. It has 546 samples in total, so what I see in it are FPKM normalised values. The first column remains the same with 60000 rows and 2 columns that is first column ENSEMBL ID and second column FPKM value.
So now the issue is I want to join all the 546 txt file into one where the first column remains the same and rest 546 column would be joined. I'm using this command:
paste -d " " *txt > merged_file.txt

This is joining the all the data sets as ENSEMBL ID and FPKM value..
But I want 1st column as the ENSEMBLE ID and rest the sample FPKM any suggestion or help would be appreciated

In R:
flist = list.files("*.txt", path=".")
fl = lapply(flist, function(x) read.delim(x, header=F))
m = sapply(fl, function(x) x$V2)
row.names(m) = fl[[1]]$V1
colnames(m) = flist

That will give you an appropriately formatted matrix named m in R, which is probably where you want to continue processing things anyway.

Let me try thanks for the quick solution  , i was trying shell script based but couldn't able

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


text:
scanBam from Rsamtools is not importing one of my reads into R
I have this read in my BAM file. It maps on chromosome 1.
I open this BAM file in IGV, and I can see the alignment on chromosome 1.
But when I open this file in R with Rsamtools:
bamContigsCel <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref     /bristolAssemblySorted.bam', param = Rsamtools::ScanBamParam(what = Rsamtools::scanBamWhat(), flag = Rsamtools::scanBamFlag(isMinusStrand = FALSE), tag = bamTags))[[1]]

I then check if the read maps to chromosome 1 in my R object but I cannot find it.
bamContigsCel$rname[bamContigsCel$qname == '000000F|arrow']
[1] II II IV
Levels: I II III IV MtDNA V X

But if I look at the BAM file, it's there. Why isn't Rsamtools importing my read into R?
$ samtools view bristolAssemblySorted.bam | grep -n '000000F' | head  -c80
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   0 0   GAAGCTGTCTAAACTTTGGC

Cross-posted on biostars

Careful, your code is using an almost certainly wrong path to the BAM file.

Note the flag; that read is mapped in a reverse-complemented manner, so isMinusStrand = FALSE is filtering it out. I tested this by making a BAM file with only that read:
@SQ SN:I    LN:1000000
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   00  *   *

Then in R:
> library(Rsamtools)
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = TRUE)))[[1]]
> blah$qname
[1] "000000F|arrow"
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = FALSE)))[[1]]
> blah$qname
character(0)

Thanks! I don't know why I was filtering out reverse-complemented alignments! Can I just use bam <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref/bristolAssemblySorted.bam')[[1]] to get all the information in the BAM file?

Yup, that should work!

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


text:
Why doesn't Biopython AlignIO.read() recognise the 'mauve' format?
On this Biopython tutorial, they describe how to import a multiple sequence alignment in the Mauve (XMFA: extensible multi fasta format). So I imported the AlignIO module:
from Bio import AlignIO
alignment = AlignIO.read(open("alignment.xmfa"), "mauve")

but I get the following error:
ValueError
Traceback (most recent call last) <ipython-input-30-396e54791d7a> in <module>()
----> 1 alignment = AlignIO.read(open("/Users/cr517/Documents/alignment.xmfa"),
"mauve")

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in read(handle, format, seq_count, alphabet)
    422     iterator = parse(handle, format, seq_count, alphabet)
    423     try:
--> 424         first = next(iterator)
    425     except StopIteration:
    426         first = None

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in parse(handle, format, seq_count, alphabet)
    366                                                 seq_count=seq_count)
    367         else:
--> 368             raise ValueError("Unknown format '%s'" % format)
    369
    370         # This imposes some overhead... wait until we drop Python 2.4 to fix it

ValueError: Unknown format 'mauve'

You're using version 1.68 or older. Mauve support was added in 1.70.


text:
Software recommendations - DNA composition
I'm looking to generate data based on the DNA composition of a region of my genomes (data is incomplete genomes from HiSeq runs in fasta format). I'm looking for software which will give me sliding windows for GC content, GC skew, codon bias etc. I would like to use something which I can use with or plug into R.

There are more R packages available that calculate GC content, for example Ape's GC.content() function.
For example:
> library(ape)
> data(woodmouse)
> GC.content(woodmouse)
[1] 0.3873347

With a sliding window is in Biostrings package.
> library(Biostrings)
> DNA <- DNAString("ACTGAAACCGTGGCAGTTTGAC")
> letterFrequencyInSlidingView(DNA, view.width=4,letters="CG")

For codon usage bias, I think cubfits or sscu packages might be useful.

Could you please [edit] your question and give an example of how the OP can use this package, what kind of results it would return, the type of input data etc? As it stands, this is more of a comment than an answer.

According to ape gc.content will only give a percentage. I'm looking for the overall distribution.

The definition of GC content is a percentage.


text:
How can I dock a protein to a nucleic acid?
I have a protein of interest and I would like to now how it interacts with RNA. I have structures of both molecules.
What tool can I use?

You mean a specific RNA molecule, right? Not any RNA in general?

Autodock Vina is a common tool for ligand binding studies, if the stretch of RNA isn't too large I suspect it may work with Vina.

I'm not sure what you meant but you can take a look at NPDock (disclaimer we wrote that tool). If you have a structure of your protein of interest, you can dock it to the structure of your DNA/RNA of interest. Mind that this is a rigid body docking which means that the structure will not change upon binding.


text:
Using data mining of papers in order to derive genomic connections
Are there any examples of researchers data mining articles and papers (e.g., from pubmed or google scholar) in order to derive relations between genes and diseases?

Do you have a specific question about text mining?

@jaslibra, it's a general question about this approach. For example to decide if a gene is tumor suppressor one. What kind of mutation are associated with pan-cancer, predict if a gene is oncogene. Basically, it has endless potential

This is really way too broad to be a good fit for this site. Unless you'd be satisfied with a yes, this can only be answered with a list fo examples. The simple truth is that everyone reads articles to learn the state of the art and then might move on to the lab to validate whatever inference they made. Please try and [edit] to make it more specific.

@terdon the question is whether a person reads it or a machine. That's the whole point.

The short answer is yes.  You could scrape pubmed for all the papers that are based on scraping data.  ;)

@0x90 both. Depends on what you need to do, that's why it's too broad. And you are not actually asking that, by the way. You are asking: Are there any examples of researchers data mining articles and papers [. . . ] in order to derive relations between genes and diseases? If you want to know how data mining works, that would be a very different question. If you want to know whether it is a machine or a human, you would have to give a specific example since both do it in different cases.

I'll provide an answer here, since the question got locked for whatever StackOverflow-ish reason. You might be interested in SmartFigures: http://smartfigures.net/ which uses data in figures to annotate an ontological database of genes, diseases, etc. See: https://vimeo.com/158555231 and https://vimeo.com/170011946 for more detail. Hope this helps and welcome to Stack Overflow!

If you feel this question was put on hold unfairly, please open a meta discussion. I am open to being corrected. I still feel this is very broad if you don't edit to limit to specific examples. I mean, of course researchers get disease associations from papers! Where else would they come from?

Also, the best way for you to reopen this question is to edit it and make it ask something specific instead of a vague "are there any examples". Explain exactly what you are thinking of (you made some more specific points in your comments, for example) do some research into what text mining is and how it can be applied and ask something specific. Or, you could change your question from "Are there any examples" to "how can this be done" which would also make it answerable.

There are many databases that have used publication scraping for oncogenic gene fusions. There are publications for the individual methods they used for their scraping and aggregation.

COSMIC - http://cancer.sanger.ac.uk/cosmic
TICdb - http://www.unav.es/genetica/TICdb/
ChimerDB - http://ercsb.ewha.ac.kr/fusiongene/ -
https://www.ncbi.nlm.nih.gov/pubmed/27899563

This existed as a closed silo, at least in 2015. Qiagen has a team of hired students and Post-Docs for collating research papers into their Knowledge Base, an extensive database that is integrated into a few of their commercial products. Qiagen's claim is that by providing a consistently-structured and well-formatted database, the process of research discovery is accelerated. I've got a few notes from a seminar they gave us a couple of years ago:

The curation is better than a free / public database because free databases typically have wrong or inconsistent data
Their knowledge base is an army of hired MDs and PhDs, 5.3M "findings" as of Nov 05 2015
Scientific articles have unstructured data, so expert knowledge is required to process the data into a "finding" for their knowledge base

I was already doubtful about their approach, and this doubt was magnified when they started getting into the details of their process. I started to lose trust and write my own thoughts about it:

The knowledge base they have created is private
The approach is to use a function (i.e. MDs/PhDs) to convert research articles into scientifically-validated associations. However, the function is effectively a black box; Qiagen allow the experts to decide for themselves how a particular article should be converted into a finding.
The findings in the knowledge base have no strength or likelihood.
The Experts were encouraged to only enter findings that were "generally accepted by the scientific community. but it seemed like only a single expert was required to get a finding into the knowledge base (i.e. findings weren't curated by other people).
When using their software (e.g. IPA), p-values of 10^-15 and 10^-22 are reported; these numbers are meaningless.
The team didn't look at methylation, phosphorylation, etc.; plans are in place to hire additional experts to look at this at a later date.

After listening to their seminar, I realised that Qiagen's idea of a general database that attempts to capture everything isn't going to work well. Experts will disagree about whether a particular association is valid, relevance (or correctness) will differ depending on the area of interest, and the established correctness of a piece of information can change over time.
It would be possible to create something that might work by adding a lot more metadata to all the "findings" in such a database, but it would need to be a huge, public, collaborative effort. Something on the scale of SRA/ENA, with encouragement from the scientific community to enter details into the database prior to publications being accepted/indexed. If the people who wrote papers entered their findings in such a database, other researchers wouldn't need to be hired to try to understand papers and guess at the findings of published research.

Thomson-Reuters is doing a similar thing, but I agree that "an association" is a complex thing and perhaps not to be taken at face value. Protein A may interact with protein B under lab conditions but whether this is a generally correct and useful thing is unclear.


text:
What kind of analysis can be done with differential expression of transcription factors?
I have two different stem cell types and their respective gene expressions from RNAseq.
I noticed that there is differential expression in some of the genes coding for histones as well as transcription factors?
What kind of interesting downstream analysis might be available?
Is it feasible to survey ENCODE for downstream TF binding sites?
Can anything interesting be investigated with regard to the histone gene data?
Thanks in advance.

TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. For histone related data, I'm not aware of how it work but perhaps you can find some leads looking at methylation studies.

Expanding my comment.
TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. Perhaps there are some databases for TF.
I found a website base on the ENCODE project with the genes regulated by each TF.You can use them for pathway or gene set enrichment analysis, but I am not aware of any specific analysis for genes regulated by a single TF.
The histone data is harder if it was methylation I would say to look if different states of the histones are in different spatial state, but with gene expression is harder to conclude anything. Nevertheless you can check methylation studies to see if gene expression is associated with methylation changes.

It sounds like you’ve already done differential expression analysis (with genes or transcripts) and a pathway enrichment analysis to find these. Differential expression of these functional groups is an interesting finding in its own right but there are several further investigations that could follow, depending on your biological research interests.
Clustering the expression patterns of the differentially expressed genes would be able to show if there are relationships between them (such as inter-gene correlations) or relationships to other phenotype data. This is commonly performed with unsupervised hierarchical clustering and displayed by a heatmap. You could also look into whether these differentially expressed genes are expressed in a tissue-specific manner. There is plenty of reference material available from ENCODE, FANTOM4/5, or GTEx that could be used for comparison (as well as cancer datasets) and the Human Cell Atlas Project (currently underway).
Another approach is pathway structure. There are many hierarchical annotation of gene functions such as gene ontology (GO terms) that overlap and are ranked: identifying smaller more specific pathways for future investigations. No pathway annotation is perfect so it’s best to replicate important findings across pathway databases (ideally different expression datasets as well). These functional groups of genes aren’t unordered: the structural relationships between pathways are well characterised in some cases. This includes many gene regulatory networks involving transcription factors. These graph structures are supported on databases such as WikiPathways and Reactome. Network analysis is relatively unexplored compared to differential expression.
While this is a bioinformatics question, there’s always a need to consider experimental validation as a future direction. If you already have a hypothesis based on your analysis, it may be more appropriate to design an experiment or contact experimental collaborators. If you have access to further genomics data, you could also explore how these differences are regulated: CAGE-Seq for enhancers and lncRNA, ATAC-Seq, and 3C/4C for chromatin structure, or Bisulfite-Seq and ChIP-Seq for epigenetic marks. An integrative analysis of expression and epigenomics data would be especially compelling and would likely account for different activity of transcription factor binding sites. Unfortunately, this can’t be done without more ‘omics data available but it would greatly benefit investigations of histones if it is feasible.


text:
Does UCSC definition for telomere making sense?
UCSC provides a database of telomere regions for each chromosome:

ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/gap.txt.gz

It reported the telomere section in chr21 as:
grep chr21 gap.txt | grep telomere
585 chr21   0   10000   1   N   10000   telomere    no
952 chr21   48119895    48129895    515 N   10000   telomere    no

Telomere is the first 10000 bases in chr21?? But ... the region are all "N" in the reference genome!
How does UCSC confidently give out telomere regions if they don't even know anything about the sequences?

Have you read any of the literature on this?

It seems so, yes.
UCSC uses a 10Kb definition and representation for telomeres in human. In the table you shared this is consistent across all chromosomes, not only chr21. This 10kb representation is an approximation to the mean value observed in newborns and is supported by evidence like the one linked below.
Taken from Wikipedia:
A telomere is a region of repetitive nucleotide sequences at each end of a chromosome, which protects the end of the chromosome from deterioration or from fusion with neighboring chromosomes. For vertebrates, the sequence of nucleotides in telomeres is AGGGTT, with the complementary DNA strand being TCCCAA, with a single-stranded TTAGGG overhang. This sequence of TTAGGG is repeated approximately 2,500 times in humans. In humans, average telomere length declines from about 11 kilobases at birth[ to less than 4 kilobases in old age, with the average rate of decline being greater in men than in women.
2,500 x 6 (TCCCAA) = 15,000
(approx. 10K)
Two interesting additional facts:
The minimum telomere length needed to ensure human telomere protective stability in white blood cells is 3.81 kb
Mouse has a 100Kb telomere representation in USCS


text:
Interpreting 0x200 flag in bwa-mem alignments
I am looking at the bwamem.h code in http://github.com/lh3/bwa and found that BWA-MEM will give flag 0x200 to what it calls MEM_F_SOFTCLIP:
#define MEM_F_PE        0x2
#define MEM_F_NOPAIRING 0x4
#define MEM_F_ALL       0x8
#define MEM_F_NO_MULTI  0x10
#define MEM_F_NO_RESCUE 0x20
#define MEM_F_REF_HDR   0x100
#define MEM_F_SOFTCLIP  0x200
#define MEM_F_SMARTPE   0x400
#define MEM_F_PRIMARY5  0x800
#define MEM_F_KEEP_SUPP_MAPQ 0x1000

But looking at the explanation in this website: https://broadinstitute.github.io/picard/explain-flags.html
It labels 0x200 as read fails platform/vendor quality checks.
If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC? Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?
The SAM/BAM specification labels 0x200 the same way as the Picard flag explaining website:

The question remains as to why BWA-MEM uses the 0x200 flag for softclipping notation rather than Fail QC.

bwa mem is using exactly the same meaning of the 0x200 flag as every other program, including picard. Don't blindly assume that that entry in the header file relates to the output flags in the SAM output. If you search the code, you'll see that bwa is using this internally as a flag for "this alignment is soft-clipped", which it handles differently for primary and secondary alignments (hard-clipping the latter).
For what it's worth, bwa mem (at least the version I just tried) doesn't set the QC failed flag. This behavour is actually pretty common among aligners, since Illumina has changed how they represent QC failure at least once.

I don’t understand this answer, as I think you’re contradicting yourself: you say that BWA-mem handles 0x200 “as every other program, including picard”, implying that it uses it for failed platform/vendor QC. But then you go on to say that it uses it for soft-clipped alignments (exactly as implied by the name in the header!). Well, which is it?

If it sets it then it has the same meaning. I just don't know off-hand if it sets it at all.

The connection between soft clipping and a read level flag is interesting.  Soft clipping occurs at a per base level and is denoted in the CIGAR string.  The failQC flag is applied to a read and usually comes directly from the sequencer. It is likely that BWA just named their internal var weird.

Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?

The failQC flag is part of the Sam Spec.

If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC?

If you expect soft clipping I would suggest trimming that sequence from the reads prior to alignment.  Also you can look at read in the input reads and see if the failQC flag is set prior to alignment to understand where it is getting set.  IGV is just displaying what is in the BAM file not making any decisions about those reads.

I wouldn't so much say that the variables are oddly named so much as it shows from the bit usage for flags (be they internal or external) that the same author wrote both bwa and SAM.


text:
Parse BAM by insertion size and get genomic coordinates
I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, file.bam.
(In terms of genomic coordinates, I just mean I would like to have a rough idea where to look at the BAM in IGV, chromosome and position.)
Based on the BAM CIGAR, I should be able to grep the BAM based on some cut-off of "large inserted" bases.
It's unclear to me how to best do this via samtools; one could parse the cigar string via samtools for insertions, but I don't understand how to translate this into the genomic position in the alignment.
This goes with various samtools wrappers, e.g. pysam in python:
import pysam

bamfile = pysam.Samfile(bamFile, "rb");

for read in bamfile:
    if not read.is_unmapped:   ## only get mapped reads
        cigarLine = read.cigar
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1: ## this means it's an insertion...
                  ## somehow filter only large insertions, and then get coordinate information

I'm not sure how to do this in a non-clumsy manner. Any help appreciated.

You basically have the right framework already in your pysam code, which I took the liberty of tidying up a bit. You just need to add a bit:
for read in bamfile:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength > 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))

Unless you have long reads, this will be close enough for IGV.

Well done! I obviously need to study the pysam API a bit more closely. Thanks for the help!

a one-liner using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html
$ java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar().getCigarElements().stream().anyMatch(C->C.getLength()>100 && (C.getOperator().equals(CigarOperator.N) || C.getOperator().equals(CigarOperator.D)))).forEach(R->println(R.getContig()+"\t"+R.getStart()+"\t"+R.getEnd()));' input.bam


text:
Quantifying reads mapping to multiple loci
I have been using STAR for our RNA-Seq samples. The final.out log file reports percentage of uniquely mapped reads along with percentage of reads that map to multiple loci (less than or equal to 10) and percentage of reads mapping to too many loci (greater than 10). However, I want to break down the multiple loci part to individual counts: Reads mapping to 2 locations, 3 locations, 4 locations .. 10 locations.
The NH tag seems to be used by STAR. However a naive read counting approach results in it reporting more number of reads than total reads.
For example, my final.out looks like this:
   Mapping speed, Million of reads per hour |       1403.36

                      Number of input reads |       53015978
                  Average input read length |       26
                                UNIQUE READS:
               Uniquely mapped reads number |       368916
                    Uniquely mapped reads % |       0.70%
                      Average mapped length |       26.45
                   Number of splices: Total |       1057
        Number of splices: Annotated (sjdb) |       0
                   Number of splices: GT/AG |       802
                   Number of splices: GC/AG |       1
                   Number of splices: AT/AC |       0
           Number of splices: Non-canonical |       254
                  Mismatch rate per base, % |       0.31%
                     Deletion rate per base |       0.00%
                    Deletion average length |       1.45
                    Insertion rate per base |       0.00%
                   Insertion average length |       1.00
                         MULTI-MAPPING READS:
    Number of reads mapped to multiple loci |       45766732
         % of reads mapped to multiple loci |       86.33%
    Number of reads mapped to too many loci |       3757890
         % of reads mapped to too many loci |       7.09%
                              UNMAPPED READS:
   % of reads unmapped: too many mismatches |       0.00%
             % of reads unmapped: too short |       5.89%
                 % of reads unmapped: other |       0.00%

Counting histogram of number of positions a read maps to using pysam:
def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts

results in
Counter({1: 330606,
     2: 86772164,
     3: 329,
     4: 38083,
     5: 31,
     6: 1094,
     7: 129,
     8: 50,
     10: 50})

The count 1 reads are fine even though they do not match the counts in final.out file since I am counting a certain category of reads (say those mapping to tRNA only), but the reads mapping to 2 locations are highly overestimated. Why is that?

You almost had the correct python code already, you just need to filter out secondary alignments:
def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        if query.is_secondary:
            continue
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts


text:
Stranded vs. unstranded library preparation protocols in RNAseq
I've been reading this paper lately:
Sailfish Enables Alignment-Free Isoform Quantification from RNA-Seq Reads Using Lightweight Algorithms
I don't really understand the second paragraph under Quantification in the Online Methods section. I quote the author:

In an unstranded protocol, sequencing reads, and hence the k-mers they contain, may originate from transcripts in either the forward or reverse direction.

As far as I understand we are sequencing cDNA here and it is a single stranded molecule so all the reads will have the same direction (unlike the DNA where the reverse complement strand is also there), so I don't see the author's point. Can someone clarify please? Thanks.

Short answer: cDNA is double-stranded.
I'll try to explain the concepts surrounding this with a nanopore read, where the cDNA component of the read (from the mouse genome) is indicated in a light grey, and the adapter components are either coloured or dark grey:

The cDNA library from this RNA was prepared by adding a polyT adapter (binding to the polyA sequence of the RNA), and carrying out reverse transcription to the start of the transcript. A strand-switching process was then used to prime and extend along the opposite strand of the DNA, creating a double-stranded full-length cDNA template. This template was then amplified using PCR to increase cDNA concentration and add in barcoding sequences.
You may notice that this particular read has a polyT head, rather than a polyA tail, attached to the VNP primer in the forward orientation. I can tell by this feature that this read comes from the reverse-complement strand of a DNA transcript, and so I know that the actual RNA transcript sequence is the reverse-complement of this read. Both the forward and reverse-complement sequences could potentially come through in the sequencer, because the adapters for sequencing are attached to both sides of the cDNA template, but the transcript direction can still be determined by looking at the orientation of the adapter sequences. When mapping this sequence to the genome, only a single strand (and its associated set of kmers) needs to be considered.
When considering what happens after mapping to the genome, there are four possible match situations:

Adapter is in reverse orientation, genome is in reverse orientation
Adapter is in forward orientation, genome is in reverse orientation
Adapter is in reverse orientation, genome is in forward orientation
Adapter is in forward orientation, genome is in forward orientation

In the case of the library preparation above, situations 2 and 3 correspond to a transcript that is encoded on the forward strand of the genome, whereas situations 1 and 4 correspond to a reverse-encoded transcript. Different library preparation methods could have a different way of encoding transcript direction. See here for more details on how these match situations can be used to produce strand-specific BAM files.
If I were to use an unstranded protocol (e.g. generating a cDNA library through other means, and attaching adapters randomly to the forward or reverse transcript strand), then the adapter direction wouldn't give me information about which strand this sequence came from. I would need to consider both strands when mapping (or more correctly, when identifying the native direction of the originating transcript), because the adapter sequence has been decoupled from the transcript direction.
On Illumina systems, this strand direction calculation is performed via the chemistry on the flow cell, and the adapters are [usually] not sequenced by the sequencer at all. Adapters from one direction bind to complementary sequences on the flow cell and are extended to produce the first read, and adapters from the other direction bind to produce the second read. Even with a stranded protocol, adapters from both directions are bound (because the cDNA template is double-stranded), but the binding order (as first or second read) can be used to determine the native direction of the originating transcript.

There are a number of library preparation methods out there. Some, like the popular dUTP-based methods, result in libraries where the orientation of a read's alignment confers information regarding the strand of the transcript from whence it came. Other protocols do more PCR cycles and don't use dUTP to allow distinguishing strands. These are then termed "unstranded" libraries and one can no longer determine the strand/direction of the originating transcript from the orientation of a read's alignment. In general, a transcript will then have reads pointing in both directions in more or less equal proportions.
You tend to not see unstranded libraries much any more. Early on they were significantly more expensive, but that's not really the case any more. There also used to be two types of stranded protocols, whereas you mostly see dUTP-based protocols these days.
BTW, you'll see these same settings in most RNAseq programs. The equivalent to an unstranded protocol setting in featureCounts is -s 0, in tophat2 it was --library-type fr-unstranded (no one should use tophat2 any more, it's just an example) and so on in other programs.


text:
Using shells other than bash
As someone who's beginning to delve into bioinformatics, I'm noticing that like biology there are industry standards here, similar to Illumina in genomics and bowtie for alignment, many people use bash as shell.
Is using a shell besides bash going to cause issues for me?

I would adjust the examples you provided. Illumina is a standard for short reads, but there are many genomics labs working mostly with PacBio or Nanopore. Bowtie is hardly a standard. Even versions 1 and 2 are very different.

@burger what do you suggest then?

No suggestion. Although I agree with all the answers so far, bioinformatics is not good with standards. Even something like a SAM/BAM file that is technically a properly defined standard that almost everyone in genomics uses has many fields that are treated differently, causing issues for a lot of tools.

A statement "this is not meant to be opinionated" doesn't help much with such a broad question like this. Do you have a particular application that you would like to use a shell for, or an indication of what "industry" you're interested in?

@burger: Do you have specific problematic SAM/BAM fields in mind?  You could raise issues at https://github.com/samtools/hts-specs/issues or at least this suggests another question to ask here…

@JohnMarshall I don't think there is a "bug" with the SAM/BAM standard. It's just that it's open-ended and different tools require different fields. I have had to modify my BAM files many times in the past because some tool expected it in a slightly different format. Technically, it's still a valid BAM before and after, but one is compatible and one isn't. If you have a BAM, you have no idea if it will work with a tool that requires a BAM file.

@burger: If you want this situation to improve, you'll have to say which particular fields you have had to modify and what the various tools' expectations were. If you do this, specs can be clarified, tools can be modified, and everyone's bioinformatics pipelines can run a little more smoothly. Otherwise, it's just FUD.

Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via #!/bin/sh or e.g. #!/bin/bash if it matters), so won't be affected by your choice of user shell.
If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell.  See Csh Programming Considered Harmful and other essays/polemics.
A Bourne-style shell is pretty much the industry standard, and if you choose a substantially different shell you'll have to translate some of the documentation of your bioinformatics tools.  It's not uncommon to have things like

Set some variables pointing at reference data and add the script to your PATH to run it:
export FOO_REF=/path/to/stuff
export PATH=/path/to/foo-x.y:$PATH
foo blah blah

These will typically be shown in Bourne-shell syntax.  By using a different shell you have to translate the export commands to your local syntax, and especially PATH munging is somewhat shell-dependent.
If you're experienced in Unix, this will be only a minor niggle.  If you're a beginner, IMHO this will add a non-negligible amount of friction on top of all the other things you're learning.

Do not use #!/bin/bash in the shebang. Having Bash installed in a nonstandard location is common enough that doing so will break often. Use #!/usr/bin/env bash instead, it should have no disadvantage.

TL;DR
SH adheres to an official industry standard, but it is not suitable for scientific computing. Bash is considered an informal standard (e.g., by Google). Bash 3 is preferable in most of situations in the world of bioinformatics.
Long answer
As already described in other answers, SH (/bin/sh, plain Bourne shell, the original UNIX shell) should fully adhere to POSIX which is a real industry standard. However, SH is too limited for scientific computing since many key features were incorporated later in SH successors, especially in Bash (/bin/bash, Bourne Again Shell): set -o pipefail, [[ ... ]], or process substitutions < () to name at least few.
In practice, it is much more difficult to write "safe" scripts in pure SH and only shell experts are usually capable to prevent unexpected behavior. For example, it may be hard to ensure that no command in a pipeline failed in the middle of computation. For Bash, various easy-to-follow defensive-programming recommendations have been developed and they should prevent all these problems. From this reason, many computer scientists, software engineers and companies use Bash as a kind of a standard. For instance, Google internal policy allows only Bash for writing any shell scripts.
Even though we cannot expect that Bash is present on completely every Unix machine (e.g., on mobile devices as @terdon pointed out), a vast majority of *nix machines used for scientific computation should have it. We also should be aware of the fact that Bash can be slower than SH and that it has recently suffered from  major security issues. Moreover, various Bash versions exist and scripts working on modern Linux machines with Bash 4 might not work on OS X, which is still based on Bash 3.
To sum up, Bash 3 is probably be the most reasonable choice for scientific computing.
Edit
I addressed the comments from @terdon and @John Marshall. In particular, I added an explanation why Bash is more suitable for scientific computing than SH (in my opinion).

Bash is not present in every Unix machine, sh is and that is not the same. Yes, Linux tends to have /bin/sh pointing to bash, but Linux is not Unix and, anyway, even in Linux /bin/sh is not always bash (Debian based systems use dash instead, for example). You can safely expect the Bourne shell (sh) to be present on a POSIX-compliant system, but not necessarily the Bourne again shell (bash).

@terdon Could you provide any reference, please? According to https://wiki.debian.org/Bash, bash is the default shell on Debian. Do you know about any (modern) *nix distro where bash would not be installed?

@terdon I will answer my question – e.g., FreeBSD. https://www.freebsd.org/doc/en/articles/linux-users/shells.html says that "Bash is not included in the default installation". Do you have an example of a Linux distribution without bash?

Some (all?) embedded Linux systems won't have bash and will instead have busybox sh. The main issue is that people tend to think that sh and bash are the same thing but they aren't. They are similar and bash is an extension of sh, but they are not the same.

@Karel: Asking about the “default shell” is ambiguous. As per https://wiki.debian.org/Shell, these days on Debian the default /bin/sh is dash while the default login shell (as listed in /etc/passwd) remains /bin/bash. This means that portable shell scripts that identify themselves with #!/bin/sh need to restrict themselves to POSIX shell facilities, while scripts that want to use bash extensions need to use #!/bin/bash.  This got tidied up the hard way a few years ago when various distributions switched to dash for /bin/sh…

@terdon @John Marshall Thank you for your comments. Compared to bash, I consider "pure" sh very limited and inappropriate for scientific computing, in particular because of some missing, but very important, features such as set -o pipefail or [[ ... ]]. My experience is that sh scripts can be very susceptible to unexpected behavior (unless the developer is a shell expert, which is usually not the case in bioinformatics). Several good and simple defensive programming strategies for scientific computing exist for bash.

This is why I would like to know whether /bin/bash might not return anything, or return a non-bash shell (I have seen such a problem only once with some obscure bioinformatics distribution).

I wouldn't do "scientific computing" in a shell no matter what shell it was. The shell should be used for, at most, handling the plumbing for basic utilities and applications. Computing should be handled by utilities and applications designed for those tasks.

@Kusalananda  How do you do scientific computing without shell? I believe that you use it at least for running your programs. If it is so, do you agree that the way how it handles errors matters?

@Karel I wouldn't do computing of any kind without a shell, but not in (with) a shell.

I’m a bit puzzled why this answer recommends the ancient Bash 3 instead of Bash 4, which itself is now almost 10 years old (announced in 2009). Bash 3 lacks crucial features such as associative arrays so it’s a severe restriction. True, macOS still ships with Bash 3 but so what? macOS is generally known to lag behind in its unix tools (and even for its Ruby and Python). Also, nitpick: it’s “Bash”, not “BASH”.

@KonradRudolph Thank you for the comment. I fixed the capitalization issue. Concerning Bash 4, I completely agree that it has many useful features. However, if it cannot be used on a substantial proportion of machines, it's a fatal problem. Whereas Python 3 can be easily installed (e.g., using Conda), upgrading Bash is complicated and easily results in serious problems. As for the associative arrays, the Google standard says the following: "If you find you need to use arrays for anything more than assignment of ${PIPESTATUS}, you should use Python."

@Karel I have a few choice word for the Google coding guidelines, none of which are acceptable in polite company. Anyway, upgrading Bash is actually trivial. Replacing it for the login shell may not be but in practice that’s unnecessary: on macOS you specify the shell in the terminal app, and other systems ship with Bash 4.

strongly agree with @Kusalananda, trying to write your pipelines entirely in shell is a mistake. There are plenty of workflow frameworks; I am partial to Nextflow, and a lot of my peers use Snakemake. Entirely shell-based pipelines quickly become unmanageable, overly complex, confusing to understand, and extremely difficult to debug. If you must use Bash, then you should aim for POSIX compliant implementations.

additionally, a lot of horrible Bash code for simpler scripts can be accomplished better with Makefiles. For beginners, you should aim to learn how to use these after you've become comfortable with basic shell scripting.

The Open Group Base Specifications Issue 7
IEEE Std 1003.1™-2008, 2016 Edition, or "The POSIX Standard" for short, is the standard that defines the interfaces and utilities provided by a Unix system.  Among these is the command line shell language and tools (see "Shell & Utilities" in the main index on the page linked above).
As far as I know, there is no shell that implement exactly what's specified by the standard, but both bash and ksh93 does a pretty good job of adhering to the standard along with their own, sometimes conflicting, extensions. The ksh93 shell in particular has had a big impact on the past development of the POSIX shell specification, but future POSIX specifications may borrow more from bash due to its wide use on Linux.
The bash shell is pretty much ubiquitous on Linux systems, and may be installed on all other Unices too. ksh93 is also available for most Unices but is usually not installed by default on Linux. ksh93 is available by default on at least macOS (as ksh) and Solaris.
If you are concerned about portability when writing a shell script (which is IMHO a good thing to be concerned about), you should make sure that you use only the POSIX utilities and their POSIX command line flags, as well as only use POSIX shell syntax. You should then ensure that you script is executed by /bin/sh which is supposed to be a shell that understands the POSIX specification.  /bin/sh is often implemented by bash running in "POSIX mode", but it may also be dash, ash or pdksh (or something else) depending on what Unix you are using.
For a Linux user, the most difficult bit in writing a portable script is often not the shell per se, but the multitude of non-standard command line flags provided by the GNU implementation of the many shell utilities.  The GNU coreutils (basic shell utilities) may, like bash, be installed on all Unices though.
Also note that bash, when running in POSIX mode (either when invoked as /bin/sh or with its --posix command line flag), is not strict about its POSIX conformance and may accept some syntax extensions to the POSIX standard.

I would not say bash as a "standard", but it is indeed likely to be the most widely used unix shell and available by default on most modern unix/linux distros. There are a few other more convenient shells like zsh that are broadly compatible with /bin/sh, but they are not as widely available. There is also C-shell and in particular its open-source implementation tcsh. C-shell is quite different from bash. Over ten years ago, I saw it was used from time to time, but nowadays, I rarely see its use, except by programmers from older generations.

The generic command sh is quite literally an industry standard, a POSIX standard, to be precise (IEEE 1003.2 and 1003.2a, available for purchase for hundreds of dollars at various websites). In theory, any script that starts with #!/bin/sh should conform to this standard. In practise, most Linux systems have a shell that is close to this standard, but has a few quirks and extensions.
Problems crop up when these quirks and extensions become standard practise in shell scripts. The Debian operating system changed to dash as their sh shell to encourage people to stop using "bashisms" in shell scripts that didn't specify a particular shell, i.e. those that began with #!/bin/sh. The dash shell tries to be as standards-compliant as possible:

dash is the standard command interpreter for the system.  The current version of dash is in the process of being changed to
       conform with the POSIX 1003.2 and 1003.2a specifications for the shell.  This version has many features which make it
       appear similar in some respects to the Korn shell, but it is not a Korn shell clone (see ksh(1)).  Only features designated
       by POSIX, plus a few Berkeley extensions, are being incorporated into this shell.  This man page is not intended to be a
       tutorial or a complete specification of the shell.

I'm not familiar with the differences, and generally try to stick to the sh man pages to instruct me with regards to correct standards-compliant shell scripts.

Note that sh is not bash. Even on systems whose /bin/sh points to bash, being invoked as sh changes bash's behavior and causes it to run in POSIX-compliant mode. The "real" sh shell (bourne shell) is something else again and not the same as bash (bourne again shell).

In Debian the default interactive shell, i.e. the one you'll be using on the command-line is bash https://wiki.debian.org/Shell yes /bin/sh will be symlinked to /bin/dash but the one people will be using live will be bash.


text:
STAR aligner multiple fastq files
I’m using STAR to align fastq files from SMART-seq2. I have raw data folder containing sub-folders with samples names the sub-folders each contain fastq file. How can I make a bash command in order for STAR to go in each of sub-folders and make alignment for the fastq file and retrieve the bam files in one folder?
Folder structure:
/plates/plate_1/rawdata/A1/A1_R1.fastq
/plates/plate_1/rawdata/B1/B1_R1.fastq
...
/plates/plate_1/rawdata/O24/O24_R1.fastq
/plates/plate_1/rawdata/P24/P24_R1.fastq

It’s 384 sample folders A1, A2, B3, F5,,, so on

Thanks I added the path structure

It seems the sequencing is single-end, and there is only one fastq per sample. If so:
staridx=/path/to/star/genome/index

mkdir /plates/star/
cd /plates/star/

files=$(find /plates/plate_1/rawdata/ -name "*_R1.fastq")

for i in ${files}
do

  sample=$(basename $i "_R1.fastq")

  STAR --genomeDir $staridx --runThreadN 8 --outFileNamePrefix ${sample}. \
    --readFilesIn $i --outSAMtype BAM

done

All fastq files will be mapped individually, and all bam files will be saved to /plates/star/.


text:
Snakemake scatter-gather with wildcard AmbiguousRuleException
My problem is when using Snakemake scatter-gather feature the documentation is basic and i modified my code according to mentioned in this link
rule fastq_fasta:
    input:rules.trimmomatic.output.out_file
    output:"data/trimmed/{sample}.fasta"
    shell:"sed -n '1~4s/^@/>/p;2~4p' {input} > {output}"

rule split:
    input:
        "data/trimmed/{sample}.fasta"
    params:
        scatter_count=config["scatter_count"],
        scatter_item = lambda wildcards: wildcards.scatteritem
    output:
        temp(scatter.split("data/trimmed/{{sample}}_{scatteritem}.fasta"))
    script:
        "scripts/split_files.py"

rule process:
    input:"data/trimmed/{sample}_{scatteritem}.fasta"
    output:"data/processed/{sample}_{scatteritem}.csv"
    script:
        "scripts/process.py"

rule gather:
    input:
        gather.split("data/processed/{{sample}}_{scatteritem}.csv")
    output:
        "data/processed/{sample}.csv"
    shell:
        "cat {input} > {output}"

I added wildcard option but, i got AmbiguousRuleException: Rules fastq_to_fasta(which is previous rule) and split are ambiguous for the file data/trimmed/Ornek_411-of-81-of-81-of-81-of-81-of-81-of-81-of-81-of-81-of-8.fasta
I tried lots of things but either rules are not calling or take AmbiguousRuleException. What am i missing, can someone help?

It's not a problem with the scatter-gather feature per-se; Snakemake is just confusing itself with what text should match what wildcards with these filenames.  You could fix this by defining a wildcard constraint for your samples, and/or changing your filename patterns to clarify what should match what.  The docs on scatter-gather show an example with different directory names and file suffixes that prevent this from happening.  Probably the simplest fix would be to follow that idea and use different directory names for the split parts, like trimmed-split and processed-split.
A hint of what's happening is the repeating ...-of-81-... thing, and the mention of ambiguity between rules.  (See also the --debug-dag option.)  If you have a sample named Ornek_411 and you want to make data/processed/Ornek_411.csv, Snakemake will look for data/processed/Ornek_411_1-of-N.csv and so on, and try to figure out how to make those.  But a filename like that could be the output of either the process rule (sample=Ornek_411 and scatteritem=1-of-N) or the gather rule itself again (sample=Ornek_411_1-of-N).  Obviously you want the former and not the latter, but it doesn't know that automatically.

Firstly thanks changing directory name pattern did trick for me. Thanks for explanation of debug-dag option too. i am aware of this logic because i got that before but with just wildcard but i couldn't set scatteritem logic


text:
Difference between computational biology, bioinformatics and biostatistics
I find that in many contexts, the terms computational biology, bioinformatics and biostatistics are often treated as functionally equivalent, and yet for students selecting PhD programs and the like the difference could be quite significant. Is there a standard or rigorous definition of these terms and the difference between them?

Biostatistics is clear, but there's a LOT of grey area between computational biology and bioinformatics.

There are many blogs and question on this on other sites, I feel the discussions always becomes rather opinion based and semantic. Here's my personal favorite blog: http://cabbagesofdoom.blogspot.ch/2015/09/is-there-important-difference-between.html

As soon as you read "biology", you can be sure there will be no accepted standard or rigorous definition ;)

I was afraid that this would descend into opinion rather than a having a straightforward answer - but I think this is a real issue that certainly students have no clear answer on, and should be addressed by this community.

NIH has had a working definition since 2000: https://www.kennedykrieger.org/sites/default/files/research_related_files/bioinformatics-def.pdf

There's also a blog post on it here


text:
How to read structural variant VCF?
The IGSR has a sample for encoding structural variants in the VCF 4.0 format.
An example from the site (the first record):
#CHROM  POS   ID  REF ALT   QUAL  FILTER  INFO  FORMAT  NA00001
1 2827693   . CCGTGGATGCGGGGACCCGCATCCCCTCTCCCTTCACAGCTGAGTGACCCACATCCCCTCTCCCCTCGCA  C . PASS  SVTYPE=DEL;END=2827680;BKPTID=Pindel_LCS_D1099159;HOMLEN=1;HOMSEQ=C;SVLEN=-66 GT:GQ 1/1:13.9

How to read it? From what I can see:

This is a deletion (SVTYPE=DEL)
The end position of the variant comes before the starting position (reverse strand?)
The reference starts from 2827693 to 2827680 (13 bases on the reverse strand)
The difference between reference and alternative is 66 bases (SVLEN=-66)

This doesn't sound right to me. For instance, I don't see where exactly the deletion starts. The SVLEN field says 66 bases deleted, but where? 2827693 to 2827680 only has 13 bases between.
Q: How to read the deletion correctly from this structural VCF record? Where is the missing 66-13=53 bases?

I'm pretty sure the END value is simply wrong in that example. Maybe send the info@ email address a note regarding this so they can fix it.

I've taken the liberty of asking the 1000 genomes project about this, I'll post the reply should it still be relevant when I receive it.

This line is wrong in several aspects. VCF4.0 is something like 5 years old and has long been deprecated. Its spec has also been moved to github. You should read the latest spec instead. The SV example there is correct.

This was fixed in the maintained VCF specifications in 2015 upon being reported as https://github.com/samtools/hts-specs/issues/84.

I just received a reply from 1000Genomes regarding this. I'll post it in its entirety below:

Looking at the example you mention, I find it difficult to come up with an
  interpretation of the information whereby the stated end seems to be correct,
  so believe that this may indeed be an error.
Since the v4.0 was created, however, new versions of VCF have been introduced,
  improving and correcting the specification. The current version is v4.3
  (http://samtools.github.io/hts-specs/). I believe the first record shown on
  page 11 provides an accurate example of this type of deletion.
I will update the web page to include this information.

So we can take this as official confirmation that we were all correct in suspecting the example was just wrong.

that was an amazing follow up! :)

So, first off, as others have pointed out, I'm pretty sure that example is just wrong. At least, the numbers don't match as you've pointed out.
That said, it is impossible to be sure without showing us the header of the VCF file as well. The INFO field (the 5th field of a VCF file) is very, very variable and depends entirely on the header lines. Each program (or human) implementing a VCF is free to choose to have whatever they feel like in the INFO field. However, each IDENTIFIER= needs to have an associated INFO line at the beginning of the file.
So, the SVTYPE, SVLEN, HOMLEN etc will have commented (start with a #) lines at the beginning of the file explaining what these values are. So check those, even though they're relatively standard, you never know, the obvious reading you used might be wrong despite its seeming so reasonable.
Here's a newer example of a VCF line for an SV taken from the current VCF specification:
##fileformat=VCFv4.1
##fileDate=20100501
##reference=1000GenomesPilot-NCBI36
##assembly=ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/sv/breakpoint_assemblies.fasta
##INFO=<ID=BKPTID,Number=.,Type=String,Description="ID of the assembled alternate allele in the assembly file">
##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">
##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=HOMLEN,Number=.,Type=Integer,Description="Length of base pair identical micro-homology at event breakpoints">
##INFO=<ID=HOMSEQ,Number=.,Type=String,Description="Sequence of base pair identical micro-homology at event breakpoints">
##INFO=<ID=SVLEN,Number=.,Type=Integer,Description="Difference in length between REF and ALT alleles">
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##ALT=<ID=DEL,Description="Deletion">
##ALT=<ID=DEL:ME:ALU,Description="Deletion of ALU element">
##ALT=<ID=DEL:ME:L1,Description="Deletion of L1 element">
##ALT=<ID=DUP,Description="Duplication">
##ALT=<ID=DUP:TANDEM,Description="Tandem Duplication">
##ALT=<ID=INS,Description="Insertion of novel sequence">
##ALT=<ID=INS:ME:ALU,Description="Insertion of ALU element">
##ALT=<ID=INS:ME:L1,Description="Insertion of L1 element">
##ALT=<ID=INV,Description="Inversion">
##ALT=<ID=CNV,Description="Copy number variable region">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype quality">
##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number genotype for imprecise events">
##FORMAT=<ID=CNQ,Number=1,Type=Float,Description="Copy number genotype quality for imprecise events">
#CHROM POS     ID        REF             ALT QUAL FILTER INFO                                               FORMAT NA00001
1      2827694 rs2376870 CGTGGATGCGGGGAC C   .    PASS   SVTYPE=DEL;END=2827708;HOMLEN=1;HOMSEQ=G;SVLEN=-14 GT:GQ  1/1:13.9

Note how the numbers do match and also note how each of the subfields in the INFO field is explained with an ##INFO line.


text:
RagTag patch error--"Tuple index out of range"
This question was also asked on GitHub
I'm trying to correct a long-read assembly with a short-read scaffold; I'm hoping to fill in the short gaps in the scaffold with the matching long-read sections. The RagTag "patch" function, if I understand correctly, is supposed to do exactly this. I use this command:
ragtag.py patch --fill-only liberibacter_central_ragtag.fasta liberibacter_central_assembly.fasta

But when I attempt to run it I get the following error:
Tue Feb  7 09:55:17 2023 --- VERSION: RagTag v2.1.0
Tue Feb  7 09:55:17 2023 --- WARNING: This is a beta version of `ragtag patch`
Tue Feb  7 09:55:17 2023 --- CMD: ragtag.py patch --fill-only liberibacter_central_assembly.fasta liberibacter_central_ragtag.fasta
Tue Feb  7 09:55:17 2023 --- WARNING: Without '-u' invoked, some component/object AGP pairs might share the same ID. Some external programs/databases don't like this. To ensure valid AGP format, use '-u'.
Tue Feb  7 09:55:17 2023 --- INFO: Running: ragtag_splitasm.py -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.agp /xdisk/baltrus/moshesteyn/liberibacter_central_assembly.fasta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Finished running : ragtag_splitasm.py -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.agp /xdisk/baltrus/moshesteyn/liberibacter_central_assembly.fasta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Running: ragtag_rename.py /xdisk/baltrus/moshesteyn/liberibacter_central_ragtag.fasta -p qseq -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.agp > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Finished running : ragtag_rename.py /xdisk/baltrus/moshesteyn/liberibacter_central_ragtag.fasta -p qseq -o /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.agp > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:17 2023 --- INFO: Writing: /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.comps.fasta
Tue Feb  7 09:55:17 2023 --- INFO: Mapping the query genome to the target genome
Tue Feb  7 09:55:17 2023 --- INFO: Running: nucmer --maxmatch -l 100 -c 500 -p /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta.log
Tue Feb  7 09:55:19 2023 --- INFO: Finished running : nucmer --maxmatch -l 100 -c 500 -p /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.ctg.fasta /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.rename.fasta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta.log
Tue Feb  7 09:55:19 2023 --- INFO: Running: ragtag_delta2paf.py /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.paf 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:19 2023 --- INFO: Finished running : ragtag_delta2paf.py /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.delta > /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.asm.paf 2> /xdisk/baltrus/moshesteyn/ragtag_output/ragtag.patch.err
Tue Feb  7 09:55:19 2023 --- INFO: Reading whole genome alignments
Tue Feb  7 09:55:19 2023 --- INFO: Filtering and merging alignments
Tue Feb  7 09:55:19 2023 --- INFO: Building a scaffold graph from the contig AGP file
Tue Feb  7 09:55:19 2023 --- INFO: Building a scaffold graph from the target/query mappings
Tue Feb  7 09:55:19 2023 --- INFO: Computing a matching solution to the scaffold graph
Traceback (most recent call last):
File "/home/u30/moshesteyn/.conda/envs/py382/bin/ragtag_patch.py", line 518, in <module>
    main()
File "/home/u30/moshesteyn/.conda/envs/py382/bin/ragtag_patch.py", line 497, in main
    match_psg = aln_psg.max_weight_matching()
File "/home/u30/moshesteyn/.conda/envs/py382/lib/python3.8/site-packages/ragtag_utilities/ScaffoldGraph.py", line 824, in max_weight_matching
    for cc in nx.connected_components(G=cover_graph):
File "/home/u30/moshesteyn/.conda/envs/py382/lib/python3.8/site-packages/networkx/classes/backends.py", line 134, in wrapper
    graph = args[0]
IndexError: tuple index out of range

Does anyone know what might be causing this error? No one else seems to have run into it on the github. I've run it with both files individually against a known reference, as well as reversing the order of the files and removing the "fill-only" option, and I have the same problem. I've looked at several of the output files, and they seem normal, at least in that they aren't empty or filled with zeroes.

Were any errors/warnings logged in any of the *.err or *.log files (in the ragtag_output folder)?

The output mentions potential AGP issues if not using -u; have you tried running it with -u included?

Can you post some assembly stats for your various assemblies? I could imagine that a filtering process under the hood removes some small artefactual contigs and exposes some sort of edge case, but it's not totally clear from what you've shown so far. For example, do you have any tiny short-read contigs that are all GGGG or something?

@Steve no errors noted in the .log, and .err is empty.

-u makes no difference, and the shortest contigs are non-repetitive between 190 and 600 bp. What kinds of stats are you looking for--spades log?

Sorry to spam the comments. Attempting to run my own patch, I've mapped the contigs individually to the reference, but for two short ones Bowtie doesn't find matches at all. Might this be causing the problem?

@schmiggle you will need to use @ to tag us when you comment, or we won't see it. Yes, any sort of stats on the assembly would be helpful, spades/quast would be fine. I doubt that the issue is a failure to find a match as that is likely a common case. I'd suggest simply filtering out contigs <500bp because they are unlikely to be very informative.


text:
How should I use dorado basecaller to calculate translocation time? Can I change models' config files?
I am trying to use dorado basecaller to process ONT data. I need to make the move table and I am doing it with --emit-move.
dorado basecaller dna_r9.4.1_e8_sup@v3.3/ /pod5_files/ --recursive --emit-moves --modified-bases 5mCG_5hmCG -v > /bams/mod_stride1.bam

The main purpose of the move table for me is to calculate the number of ONT device measurements before basecaller emits one base (calculate translocation time per base). But the stride value in config.toml file is fixed to 5 (stride=5). I need to calculate the instantaneous translocation time per base (in seconds), because I am looking to find a measure similar to IPD values in PacBio, when we can get the time for sequencing each base.
I was wondering to change the stride value in config.toml file in the model folder, from stride=5 to stride=1, to directly infer the number of measurements without using Tombo for correcting and re-squiggling.
Unfortunately, even though this command works perfectly fine when the stride=5, when I set stride=1 it always gives me the following error: "CUDA out of memory", even when I triple the size of RAM for the job.
My questions are:

Does setting stride=1 make sense for this purpose? Or do I have to do something else? (I know we have some K-mer models too)
Is this even possible to change the stride value on config.toml file?
How should I figure out how much RAM I need for the job?

In pod5 files the number of samples prior to the adapter/sequence should be encoded in the metadata, and can be looked at using ONT's pod5 inspect tool:
$ pod5 inspect reads pod5_file.pod5

read_id,channel,well,pore_type,read_number,start_sample,end_reason,median_before,calibration_offset,calibration_scale,sample_count,byte_count,signal_compression_ratio
00445e58-3c58-4050-bacf-3411bb716cc3,908,1,not_set,100776,374223800,signal_positive,205.3,-240.0,0.1,65582,58623,0.447
00520473-4d3d-486b-86b5-f031c59f6591,220,1,not_set,7936,16135986,signal_positive,192.0,-233.0,0.1,167769,146495,0.437
...

The average translocation speed can be worked out by comparing the sample_count to the number of called bases.
More details here.
Unfortunately, due to the way the basecaller models work, I don't think there's any easy way to create a 1-sample lookup table for precise instantaneous translocation calculations.
If it is just an out of memory error, then it should work (but very slowly) if you set the basecaller batch size to 1 (i.e. -b 1). If even that doesn't work, something else is going on and it's unlikely to be a memory problem.


text:
Problem of writing scilab code on a project - Morris-Lecar model
This is a computational neuroscience project. We are asked to plot the Morris-Lescar model (question 3 below).
The model should be plotted with scilab or R.
however, I don't know where to start? For the moment, the problem is more in writing the code to represent the equations and systems of equations than in actually solving them.
Any clues as to how to start writing our code would be a huge start!

In concrete terms, we don't know where to start with my working group, even for question 3.a..
We know that nullclines are associated with $dV/dt = 0$ and $dW/dt = 0$, which gives the following entries:
$W = \frac{I-g{Ca}M{\infty}(V)(V-V_{Ca})-g_L(V-V_L)}{g_K(V-VK)}$
$W = W{\infty}(V)$
But here for example, no idea how to write it on scilab, or how to display them graphically by varying I.
I hope the application is now clearer.
Appendix: scilab code.
E_L = -50;
E_Ca = 100;
E_K = -70;
V1 = -1.2;
V2 = 18;
V3 = 2;
V4 = 30;
g_L = 2;
g_Ca = 4.4;
g_K = 8.0;
C = 20;
gamma = 0.04;

minf(v)=.5(1+tanh((v-va)/vb))

winf(v)=.5(1+tanh((v-vc)/vd))

tauw(v)=1/cosh((v-vc)/(2*vd))

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

I tried to make the question more to the point, but I think would be easier to answer if you would be even more specific.


text:
Error occurrence after merging files with bcftools: wrong number of fields?
This question was also asked on Biostars
I have multiple vcf of CASES and CONTROLS variations annotated by VEP, SNPEff and SnpSift.

First pair vcf -> only variations| CASES and CONTROLS
Second pair vcf -> variations + SnpEff | CASES and CONTROLS
Third pair vcf-> variations + SnpEff+VEP+SnpSIFT| CASES and CONTROLS | at so on
The data has non-overlapping samples. I tried to merge corresponding vcf of CASE and CONTROLS into one file with bcftools merge.

The code is:
 bgzip -c file.vcf > file.vcf.gz
 tabix -p vcf file.vcf.gz
 bcftools merge -o merged.vcf.gz -Oz file.vcf.gz file_2.vcf.gz

Unfortunatelly in files with SnpSift annotations i got this error:

Error at chr1:20724637: wrong number of fields in dbNSFP_ada_score?

Please anybody knows how to solved this error or what does it mean?
grep header:
INFO=<ID=dbNSFP_ada_score,Number=A,Type=String,Description="Field 'ada_score' from dbNSFP">

-G on chr1:
CHROM POS ID REF ALT QUAL FILTER INFO
chr1 20724637 rs1280663306;rs1445195248 TGGGAGGGAGGGAGAGAGGT TGGGAGGGAGAGAGGT,TGGGAGGGAGGGAGGGAGGT,TGGGAGGGAGGGAGAGAGGG,TGGGGGGGAGGGAGAGAGGT,TGGGTGGGAGGGAGAGAGGT,TGGGAGAGAGGGAGAGAGGT 2.42927 . ANN=TGGGGGGGAGGGAGAGAGGT|splice_region_variant&intron_variant|LOW|SH2D5|SH2D5...SNV|1||||||||||||||||||||||||||||||||||||||||||||;dbNSFP_ada_score=1.43020658327008E-4,0.00157457299236331,5.16496450143413E-5;dbNSFP_rf_score=0.016,0.054,0.0;dbNSFP_PHRED=5.632,0.077,6.365,7.459,2.716
chr1 20724637 rs1280663306;rs1445195248 TGGGAGGGAGGGAGAGAGGT TGGGAGGGAGAGAGGT,TGGGAGGGAGGGAGGGAGGT,TGGGAGGGAGGGAGAGAGGG,TGGGGGGGAGGGAGAGAGGT,TGGGTGGGAGGGAGAGAGGT,TGGGAGAGAGGGAGAGAGGT 2.42927 . ANN=TGGGGGGGAGGGAGAGAGGT|splice_region_variant&intron_variant|LOW|SH2D5|SH2D5|transcript|NM_001103161.2|protein_coding|5/9|...||||||||||||||||||||||||||;dbNSFP_ada_score=1.43020658327008E-4,0.00157457299236331,5.16496450143413E-5;dbNSFP_rf_score=0.016,0.054,0.0;dbNSFP_PHRED=5.632,0.077,6.365,7.459,2.716

dbNSFP_ada_score is an ensemble prediction score, as calculated by dbNSFP.


text:
How can I add a ";" at the end of each fasta header line in a file?
I have a large fasta file with RNA sequences. I need to add a ; at the end of each header line (lines beginning with >) in the file.
For example, I have:
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

And, my desired output is:
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

One way, using sed and a regex address to select only the header lines and apply a substitution to append a semicolon to the end of each line:
sed '/^>/s/$/;/' in.fa > out.fa

Or, to edit the file in place:
sed -i '/^>/s/$/;/' in.fa

Or one way using awk and a ternary construct to print the line followed by a semicolon or an empty string, depending on whether or not the line is a header line:
awk '{ print $0 (/^>/ ? ";" : "") }' in.fa > out.fa

Alternatively, you may prefer:
awk '/^>/ { print $0 ";"; next } { print $0 }' in.fa > out.fa

Which is the same as:
awk '/^>/ { print $0 ";"; next }1' in.fa > out.fa

If you have perl installed, all you need is:
perl -pe 's/$/;/ if /^>/' in.fa > out.fa

... using the legendary Perl pie
perl -p -i -e 's/^(>.*)/$1;/' mybacteria.fa

Input
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG
A

Output
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

Details I copy pasted your first code into a file mybacteria.fa then ran the Perl pie script. Perl pie is very quick, there alternative regexs (regular expressions) that would be slightly more efficient, but Perl pie is fast regardless. In fact you don't need to use regex in this example and would be faster via a if $line == '>'; declaration (certainly true for Pythonic equivalents if line == '>') but it doesn't offer the ease and flexibility, which is combined with speed, that Perl pie does. The code will work in Linux (any flavour/flavor) or Mac OSX Darwin (i.e. Terminal app)

It will work on any flavour/flavor so long as Perl is installed. Note that not all distributions ship with Perl pre-installed. Also, the capture group is not necessary - all that's needed here is: perl -pe 's/$/;/ if /^>/' in.fa > out.fa. Note that using -i without an extension will modify the input file in-place without making a backup of the original.

In fairness I did point that out in the discussion below the code. The efficiencies are generally not noticeable. However, most regex manipulations take place within the string rather than at the end requested here, and capture can be one way to deal with different situations beyond this one.

Fair point. Perl is included in most popular Linux versions including Ubuntu, Debian and Suse. Perl isn't preinstalled on RedHat or freeBSD after release 5 and nobody uses Kali in bioinformatics. Let me know if there's anything missed.

Thank you very much -- this worked well and was simple to execute.

@CallahanMcGovern if an answer worked for you, can you please mark it as accepted?

seqkit replace -p "$" -r ";" my.fasta > my.renamed.fasta

Where -p is the pattern to match the end of the header text $
Where -r is the replacement pattern

https://bioinf.shenwei.me/seqkit/usage/#replace
Seqkit is an extremely powerful toolkit for working with fasta/fastq files. It can a lot with a single cli command.

I recommend to use Biopython as a clearer way.
from Bio import SeqIO

records = []
for record in SeqIO.parse("file.fas", "fasta"):
    record.description += ";"
    records.append(record)
SeqIO.write(records, "file.fas", "fasta")

Before:
>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA

After:
>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA

Clearer than what? Why is using a complex package requiring several lines of code "clearer" than using standard *nix tools and one liners?

@terdon For the same reasons fasta, nexus, &c. are explicit. It is easier to read at a glance.

I find this much harder to read than the standard tools with their simple syntax. I need to grok 6 lines of code here and manually enter the file name twice with no autocompletion. You also don't mention that you need to install biopython, while the other solutions work out of the box on pretty much any *nix system. On the whole, this is longer, more complex, less portable and significantly slower than either the perl or sed one-liners. This took 6.7s compared to 2s for Perl and 1.6s for sed on a ~300MB file. It does work, but let's stick to the technical and remove the opinion :).

@terdon Of course, you are right from this side. But it reminds me the times when I have only 105 bytes of programme memory & try use the address byte of GOTO as hex-command forwarding on it from another GOTO because it was more effective :-)

Heh, fair. But I just honestly don't find this in any way clearer. It certainly took me more seconds to understand this than any of the other answers. You are obviously more used to python than to perl or sed, so of course you will find it clearer. I am familiar with python, perl and sed, and of the three solutions, this was the least clear since it had more text to read. Everyone's experience is different! :)

@terdon As someone who has taught programming to juniors I would bet a significant amount of money that a larger percentage of 'people in IT' will correctly "decode" what this does compared to any regex based solution. The strength of regexes is their flexibility, not their cleanness or readability. Appending a character at the end of a line really isn't playing to the strengths of regexes.

@DavidMulder oh sure, programmers are more likely to be familiar with python, no question about that. Old school bioinformaticians/*nix geeks like myself, less so. In any case, it was the unsubstantiated "clearer" that rubbed me the wrong way, since it is so very subjective and happens to be wrong in my own case.  And it also feels like using a tank to swat a fly: we're loading all sorts of fancy tools that just make this a slower, less efficient solution.

@DavidMulder I'm with @terdon on this one - I don't find this solution clearer at all. It might appear easier to read at a glance, perhaps because record.description += ";" seems rather innocuous, but what happens if we're now presented with a FASTA file without a description field? Well, you need to know that in this case the description attribute is actually assigned the record id/name. So appending a semicolon to the description attribute will actually create a description field in the output file, basically duplicating the record id/name. Probably not what anyone expected or wanted.

Steve's point deserves emphasis. I love Biopython but its implementation of the id versus description attributes on SeqRecord objects never fails to trip me up.

This script loads all sequences into memory, which wouldn't work for huge input.


text:
Introducing a fixed number of random substitutions in a sequence
I'm writing a function that introduces $n$ substitutions inside a sequence of nucleotides.
I have a working version, but I'm looking for any other ways of doing this:
from random import choice as random_choice
NUCLEOTIDES = set(["A", "C", "G", "T"])
def randomly_substitute(seq: str, num_mutations: int = 1) -> str:
    if num_mutations > len(seq):
        raise ValueError(
            f"{num_mutations} is greater than the length of the input sequence"
        )
    index_choices = list(range(len(seq)))
    mutated_seq = seq
    while num_mutations != 0:
        index_choice = random_choice(index_choices)
        index_choices = [el for el in index_choices if el != index_choice]
        chosen_nucleotide = seq[index_choice]
        if chosen_nucleotide not in NUCLEOTIDES:
            raise ValueError("Not a nucleotide: {chosen_nucleotide}")
        possible_mutations = [el for el in NUCLEOTIDES if el != chosen_nucleotide]
        chosen_mutation = random_choice(possible_mutations)
        mutated_seq = (
            mutated_seq[:index_choice]
            + chosen_mutation
            + mutated_seq[index_choice + 1 :]
        )
        num_mutations -= 1
    return mutated_seq

In particular, the function maintains a list called index_choices that is of size len(seq), and keeps re-building this list after each introduced substitution. This feels not very efficient.

Here's a slightly different approach that should be more efficient. The core of the idea is to have two functions, one which changes a given nucleotide to something other than what it was, and another to perform the substitution on the input sequence. This way, you only need to build the sequence once.
#!/usr/bin/env python

import random

def change_to_not_self(old, nts):
    if old not in nts:
        raise ValueError(
            f"{old} is not a valid nucleotide"
        )
    ## Pick the new nucleotide from the list of acceptable
    ## values minus whatever letter we have passed as 'old'.
    new = random.choice(list(set(nts) - set(old)))
    return new

def randomly_substitute(seq: str, num_mutations: int) -> str:
    if num_mutations > len(seq):
        raise ValueError(
            f"{num_mutations} is greater than the length of the input sequence"
        )
    NUCLEOTIDES = set(["A", "C", "G", "T"])
    seq = list(seq)
    indices = random.sample(list(range(len(seq))), num_mutations)

    for index in indices:
        seq[index] = change_to_not_self(seq[index], NUCLEOTIDES)
    return ''.join(seq)

sequence = 'ACTGACGGTACGTACAACG'
mutated_seq = randomly_substitute(sequence, 5);
print(mutated_seq)

Nice!
Although when you say There is nothing there to check if index_choice has been used in a previous iteration of the while loop,, that's not correct IMO, as that's what line index_choices = [el for el in index_choices if el != index_choice] does. Your solution is clearly an improvement though, so approving it :))

@bricoletc thanks :) My point was that you were setting index_choice = random_choice(index_choices) and only then setting index_choices = [el for el in index_choices if el != index_choice], so in the unlikely case where 3 iterations of the loop all chose the same index in random_choice(index_choices) you could have changed the same position 3 times. So chosen_nucleotide could be repeated.

How? We start with all possible indices, say [0, 1, 2, 3]. Then we iterate num_mutations time a process where we select one index, say 0, remove it, so index_choices is [1,2,3], then we select another index, say 2, and index_choices becomes [1,3], etc. No?

Ah, duh, of course! You're absolutely right, @bricoletc, my bad.

In particular, the function maintains a list called index_choices that is of size len(seq), and keeps re-building this list after each introduced substitution. This feels not very efficient.

Yes I agree with this. Personally, I'd simply use,
 index_choices = random.randint(1, len(seq))

Lists are expensive in Python and you're creating a list for every value of the sequence when randint is quicker and random.choice is if the numbers are not sequential (there's absent values). Here you just want a rectangular distribution from a straight numerical range.
mutated_seq is perfectly good. A numpy solution is would speed this up BTW. I can see where you've got the idea of random.choice from, but its not needed here.

But, index_choices is a list of the remaining possible choices; if num_mutations > 1, then I can't choose any int between 1 and len(seq), because that might choose the same position twice

Strings in Python are immutable (can't be changed), so if you want to do something like introducing random mutations it may be better to use an intermediate mutable object like a list. Here's one way to do it via sample and choice, which will raise an error if the number of samples desired is greater than the sequence length:
#!/usr/bin/env python
from random import sample, choice
OTHER_NUCLEOTIDES = {"A": ["C","G","T"],
                     "C": ["A","G","T"],
                     "G": ["A","C","T"],
                     "T": ["A","C","G"]}

def randomly_substitute(seq: str, num_mutations: int = 1) -> str:
    seq = list(seq)
    for pos in sample(range(len(a)), num_mutations):
        seq[pos] = choice(OTHER_NUCLEOTIDES[seq[pos]])
    return ''.join(seq)
```