text:
How do I rename fasta headers?
I am trying to use the seqkit replace command to replace chromosome names in the format chr_I, chr_II, ... to I, II, .... I am using the following command:
seqkit replace -p "(.)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

My keyValues.txt file contains the following:
chr_I   I
chr_II  II
...

The two columns are separated by a tabulation.
I get the following output:
[INFO] read key-value file: keyValues.txt
[INFO] 6 pairs of key-value loaded
[ERRO] pattern "(.)" matches multiple targets in "chr_I", this will cause chaos

I don't know what I am doing wrong. I chose the (.) pattern to match the whole header but it seems to be wrong. Any help would be appreciated.
Update
All the headers are shown below:
grep '>' mySequence.fasta
>chr_I
>chr_II
>chr_III
>chr_IV
>chr_V
>chr_X

Show some of the headers, the pattern argument is probably not suitable. I am guessing that you'll get headers with duplicate names now (which will cause chaos).

I have written all the header names in my question.

There are 4 matches to your pattern "chr_I"

Oh, so the key is a regular expression? I tried putting a dollar at the end of the key names, but I get the same error.

Why not use sed instead?
sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta

Or even simpler:
sed 's/chr_//' mySequence.fasta > mySeq.fasta

Worked. I just don't know sed. I need to learn it.

Great! It works a bit like grep (which you know how to use), but then with replacement function.

I would suggest using sed 's/^chr_//' instead. Anchoring the pattern to the beginning of the line should speed things up for large files. Otherwise, sed will have to scan the entire length of each line in case there's a match.

@terdon, good point. But does it work with fasta? Remember the > at the beginning of the header.

@b.nota Duh! I should learn to read :) Indeed, sed 's/^>chr_/>/ would be better.

Minor correction for the last: sed 's/^>chr_/>/' ( in order to keep the >)

I like using awk when there's any kind of lookup involved: something like
awk 'FNR==NR {
    hash[">"$1]=$2;
    next;
}

if ($1 ~ /^>/) {
    print ">"hash[$1];
} else {
    print $1;
}' keyValues.txt mySequence.fasta

Basically the FNR==NR check tells awk to work only on the first file (useful primer here), keyValues.txt and create an association for each key (eg. chr_I) with its value (eg. I). The rest of the code after the next works only on mySequence.fasta, printing out the lookup value only if the line is a fasta header, as checked by the $1 ~ /^>/ condition.

Here is an alternative to sed: cat test | tr -d "chr_"
The pipe is to output the file, while the tr command is to translate and change the file

The use of the cat command is completely redundant here (and imposes a runtime overhead). Just do tr -d 'chr_' < test.

Your issue is that the pattern . matches a single character. To match the whole header, use .*:
seqkit replace -p "(.*)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

Alternatively as pointed out by others you may not really need keyValues.txt here as you're just removing the chr_ prefix:
seqkit replace -p "chr_(.*)" --replacement '$1' mySequence.fasta

(Don't forget the single quotes for $1 else the shell will expand it)


text:
Correct way to assemble reference transcriptome - what's --samples_file option in Trinity?
I have 101 samples from 9 tissues of a nonmodel species, coming either from SRA or sequenced by my lab. I want to generate a reference transcriptome with Trinity for further analysis. In order to do this, I downsampled the reads in order to have 20 to 30M reads for each tissue, for a total of 270M PE reads, which I combined in two files, one for left reads and one for right ones.
Now,  I came across the --samples_file flag in Trinity and I'm not sure if I should use it instead of concatenating the transcriptomes. Any help on good practice on this?

--samples_file flag is for DE analysis, you can put info about biological replicates in that file. You say that you want to make a ref transcriptome, so you don't need a samples file for DE analysis.
Why do you down sample? Are you not afraid you will not have enough reads to have a complete transcriptome now?

Thanks, it's what I was thinking. // I have 270 million reads, I think it's enough. Using all transcriptomes means having >750 million reads, I don't have the computational resources to assemble it.

Okay, I understood 270 PE reads...

Yes, I actually wrote 270 PE reads, my bad.

source: Stack Exchange

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

sort: votes


text:
Why do NEBNext indexing primers have sequence between the p5 oligo and index?
In a previous post I asked Why do NEB adapters have non-complementary sequence?
Since then, I realized that there is some other sequence in the p5 indexing primer, as well as in the p7 indexing primer.
Here is a diagram of the NEBNext protocol. The parts that I am confused about are in step 5.

My first question is about the p5 index. Why is there a sequence 5'GATCTACAC 3' between the p5 flowcell oligo and the index?

My second question is about both indices. Why is there no gap between the p7 index and the flowcell annealing oligo?

For reference, I fetched these indexing primer sequences from NEB on page 18 and 19 of NEBNext Multiplex Oligos for Illumina instruction manual.

Question 1: The additional sequences is needed because that is complimentary to the P5 sequence anchored to the flow cell. It is also the site for priming the Index 2 read on a MiSeq. The additional sequence is also needed for the read 1 primer in the cartridge to anneal to the correct place for the molecules.
Question 2: The "GT" sequence you're referring to is necessary for the Index 1 and Read 2 primer to anneal to the correct place.

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


text:
What are some good sources of 3rd order gene expression data
I need good sources of gene expression data in the form of 3rd order tensors. Typically the commonly available datasets are in the form of a matrix, for instance, $sample \times gene$ or $gene \times time$. I need good sources for 3-dimensional data, for example $gene \times sample \times tissue$ or $gene \times sample \times time$ ... etc.

Analysis packages rarely use that format for analysis, so find a nice time course dataset (there won't be many matching your needs) and reshape the matrix for it into the correct dimensions. As a aside, it's rarely useful to refer to a 3D matrix or 3 factor experiment as a "3rd order tensor" outside of machine learning.

I am working on tensor decomposition. I am interested in testing an idea on gene expression 3-d dataset. The papers who have worked on such gene expression dataset haven't made their data public. So...

@Satwik, you can ask them for the data and do a collaboration (if they and you wish).

I familiar with the inflammatory bowel disease  which is quite complicated:
It affects multiple sites and depending on when it appears it is different. So you can look at the GEO for microarrays and RNA-seq data of this disease and I'm sure you'll find datasets from the same patients from several tissues or time.
In this disease there is a multi'omics project that fits your needs. You can find (and download) data of high dimensionality: there are samples for the same time, tissue (but different locations), gene, over a period of time and for several omics (not only expression data).
The only problem is that there are some mistakes on the metadata.
However, I don't think they are usually described as 3rd order tensors and I doubt you'll find enough data to train a machine learning method.


text:
How to calculate logCPM across all samples?
Using edgeR for differential analysis between Tumor and Normal gave me differential expressed genes with logFC, logCPM, PValue and FDR.
From the details of glmTreat function I see that logCPM is average log2-counts per million, the average taken over all libraries.
The output table I got from glmTreat:
Geneid     logFC  unshrunk.logFC    logCPM  PValue       FDR
Gene1   4.985298235 5.573931983 1.236660486 1.55E-54    1.87E-50
Gene2   4.613438634 5.52126484  0.920343233 2.59E-53    1.56E-49
Gene3   5.250601296 5.356432653 2.751294034 2.48E-50    9.95E-47
Gene4   4.943049159 5.379741165 1.361393757 2.06E-45    6.19E-42
Gene5   6.100121754 6.117580436 5.774385315 2.95E-43    7.11E-40
Gene6   4.722697891 5.320461275 1.120685402 3.05E-42    6.11E-39
Gene7   5.246129853 5.497012001 1.902992053 2.78E-40    4.78E-37
Gene8   3.878773277 4.956636276 0.776208741 8.66E-39    1.30E-35
Gene9   4.441752496 4.94930499  1.132652682 4.37E-38    5.83E-35

To calculate logCPM given in the above table manually I did like this:
logCPM <- cpm(y, prior.count=2,  log=TRUE)

I got the logCPM values for all the genes. Then I calculated Average across all the samples for each gene. I got like below:
Geneid  Average
Gene1   0.686560246
Gene2   0.617115826
Gene3   1.075975225
Gene4   0.692050878
Gene5   1.277556065
Gene6   0.638358189
Gene7   0.689323163
Gene8   0.60700396
Gene9   0.662115092

Why there is so difference? Where I'm doing wrong?

Did you also try aveLogCPM(y, normalized.lib.sizes=TRUE, prior.count=2, dispersion=NULL, ...)?

glmTreat() literally uses aveLogCPM under the hood in fact.

@DevonRyan Yes you are right. thanks. It is calculated across all samples. I want to calculate aveLogCPM for each condition differently (Tumor and Normal separately) How to do that? Just taking the raw counts of Tumor condition and using aveLogCPM function on it is fine?

If you want the values by group then subset y to contain the samples of interest and then feed that to aveLogCPM().

yes I took the matrix of only Tumor samples and genes with raw counts and applied aveLogCPM()

This post on the Bioconductor mailing list explains the difference: https://support.bioconductor.org/p/100463/

Thanks for your response. Could kindly provide a summary of the link here? Link answers are generally not encouraged.

In short, Aaron Lun explains that the output of aveLogCPM is closer to
log2(rowMeans(cpm(y)) (where the log-transformation is done after the averaging), than to:
rowMeans(cpm(y, log=TRUE))


text:
Which data is being used for violin plot?
Sorry I just got totally confused conceptually.
If this is my raw count data
Likely, Seurat divides each value by sum of the column afterward times by 10000.
which gives so
                   s1.1        s1.2 s1.3       s1.4 s1.5
DDB_G0267178 0.00000000 0.009263254    0 0.01286397    0
DDB_G0267180 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267182 0.00000000 0.000000000    0 0.03810585    0
DDB_G0267184 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267188 0.02640801 0.000000000    0 0.01286397    0

Literally Seurat claims that final step would be taking natural log of the above matrix as I tested with
library(seurat)
mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
mat

    [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    1    1    1
[2,]    1    0    2    2    3
[3,]    1    1    1    1    2
[4,]    1    1    0    1    0
[5,]    0    0    0    1    3
mat_norm <- LogNormalize(data = mat)
mat_norm
[1,] 7.824446 .        7.824446 7.419181 7.014015
[2,] 7.824446 .        8.517393 8.112028 8.112028
[3,] 7.824446 8.517393 7.824446 7.419181 7.706713
[4,] 7.824446 8.517393 .        7.419181 .
[5,] .        .        .        7.419181 8.112028

My confusion is: seurat@data dose not give me natural log transformed data rather returns only data divided by column sum and scalled by 10000.
So when plotting a violin plot if data are log normalised plot should not look as below picture
When I used
VlnPlot(object = y, features.plot = "DDB_G0277853", x.lab.rot = FALSE, y.log = TRUE)

plot changed totally. Could someone please help me in getting idea about seurat@data that is not really log transformed and by which data violin plot is being produced??? So if a gene has 10000 read counts, without log transformation violin plot should be so
This is my seurat object

The normalized and log-transformed values are used for the violin plot.
The argument y.log changes only the display of the data (scaling of the y axis).
Seurat has very good documentation. Section 7 in the FAQ explains what data is stored in the object:

How is data stored within the Seurat object? What is the difference between raw.data, data, and scale.data?

raw.data The raw data slot
(object@raw.data) represents the original expression matrix, input
when creating the Seurat object, and prior to any preprocessing by
Seurat.
data The data slot
(object@data) stores normalized and log-transformed single cell
expression. This maintains the relative abundance levels of all genes,
and contains only zeros or positive values. See ?NormalizeData for
more information. This data is used for visualizations, such as violin
and feature plots, most differential expression tests, finding
high-variance genes, and as input to ScaleData (see below).
scale.data
The scale.data slot (object@scale.data) represents a cell’s relative
expression of each gene, in comparison to all other cells. Therefore
this matrix contains both positive and negative values. See ?ScaleData
for more information.

Thank you, you are definitely alright. The problem is; when I am normalising values for some genes and compare them with output of seurat@data, some genes are ok, but some genes seems to be without log normalisation and all stored in seurat@data. I guess maybe after calling by 10000, only some especial range of expression values go through the log normalisation, for example seurat@data gives me 0.0026 for a gene in a cell, when I calculated  by pen, I noticed this value is just deciding the raw count on column sum multiply by 10000. Please correct me if I am wrong

In your first example, the last row does not match: DDB_G0267186 vs DDB_G0267188. Note that object@raw.data contains all cells, meanwhile object@data contains only the filtered cells. That might explain the difference.

@FereshTeh look for comparison at the formula Seurat uses: log1p(value/colSums[cell-idx] *scale_factor) (note the implicit +1 and logn in the formula)

Those two graphs don't look radically different.   The y axis is definitely not log scaled on the top one, and it looks like it is scaled in the second one (y.log = TRUE)

Sorry, what you think about logNormalize function in Seurat that does not return log  transformed data?

Sorry, for instance, for gene 1 in cell 1 I have 1040 RAW read counts. the sum of raw read counts in cell 1 for all genes is 373762. Based on the seurat documentation for normalisation, I should do so;  ln ( 1040/373762 * 10000) that gives me 3.6 that is right. So in top violin plot for this gene seurat has used the natural log transformed data as stated and in seconds violin plot I am just taking another log of y axis. But question is: why seurat@data does not return values as I expect and I called for one gene rather give me only scaled data by 10000. could you help me?


text:
Human Genome Alignment and Variant Calling Benchmarks
I am doing a project that requires bench-marking the performance of computational tools for alignment and variant calling of human genome sequencing data. In particular I would like to establish benchmarks comparing the performance of GATK and ADAM/Avocado.
My main question is which genomic data-sets should I use as the input to the data processing pipeline. In Frank Nothaft's (author of ADAM and Avodaco) thesis he reports bench-marking the processing of NA12878 from the 1,000 genome's project in Avocado as compared to GATK's HalotypeCaller. He also reports benchmarking on 270 samples from the Simons Genome Diversity Data-set.
I am interested in reproducing this portion of his work, but do not know where I can get a hold of these datasets. Where can I get a copy of the NA12878 genome (BAM file) and also were can I get aligned samples from Simons Genome Diversity Data-set?

NA12878

My recommendation: download raw reads from Illumina BaseSpace and do alignment yourself. Aligning ~40X worth of human reads takes overnight with ~16 CPU cores.
You can acquire Platinum Genome data from ENA. There is a BAM available for download, but that was done by Illumina's own aligner against a wrong reference genome. You'd better do alignment by yourself. Another problem with this dataset is it is too good to be representative. Modern data are worse.
GIAB provides BAMs here:
ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/

It is not "typical", either.

SGDP

Instructions on downloading SGDP BAMs are available on this page, in particular this Word document. Warning: you need ~30TB disk space for these BAMs. I also heard it is tricky to download data.
There is apparently a mirror from cancer genomics cloud. I don't know how easy it is to access.

I'm having some trouble finding the fastq files for NA12878 on Illumina BaseSpace. Do you have any suggestions on where to find these?


text:
Ordering the cells by the expression of specific genes
I have 209 cells. I want to show the expression of two genes by ordering my cells. Imagine something like a plot on which y axis shows the read counts range from 0 to 10000 and x axis shows the number of cells (there is a link of an example picture that I wish to produce the same). So that two colours bars each colour for a gene.
I have used this code but result is not as I want at all
# Library
 library(tidyverse)

 # Create data
 value1=abs(rnorm(26))*2
 data=data.frame(x=LETTERS[1:26], value1=value1, value2=value1+1+rnorm(26, sd=1) )

 # Reorder data using average?
 data = data %>% rowwise() %>% mutate( mymean = mean(c(value1,value2) )) %>% arrange(mymean) %>% mutate(x=factor(x, x))

 # plot
 ggplot(data) +
   geom_segment( aes(x=x, xend=x, y=value1, yend=value2), color="grey") +
   geom_point( aes(x=x, y=value1), color=rgb(0.2,0.7,0.1,0.5), size=3 ) +
   geom_point( aes(x=x, y=value2), color=rgb(0.7,0.2,0.1,0.5), size=3 ) +
   coord_flip()

This is link of my data and my result picture

And to some data

Could you add an example of what a lollipop is? (a link or something like that?) And what have you tried to make that plot (To avoid suggesting to do the same and help you more specifically)

Thank you, I edited my post. Actually lollipop plot itself dose not matter that can be replaced by any other plots. The problem is that I could not figure out how to plot the same as picture and what I have tried up to now did not answer my question :(

could you further explain how this picture doesn't help you ?. I am still trying to figure out how I can help

Thank you, because on gene seems to be compressed at the end  and does not allow to another gene to be cleared. I don't know how to have another gene on another axis or showing y axis in log10

So, your problem is how to convert an axis to log10? This is answered in Stack Overflow and other venues. In ggplot2 you can't have two axis by design: it is against the grammar of graphics.


text:
How to superimpose barplots in karyotype graphic using the chromPlot Bioconductor package?
I have a data frame containing 4 columns:

Chrom: chromosome name of the feature
Start: starting base of the feature
End: end base of the feature
Name: name of the feature, e.g. 5S ribosomal RNA

The head of the data frame is shown below:
      Chrom    Start      End                        Name
    1     1 15390242 15390720 16S ribosomal RNA (partial)
    2     1 15390900 15391794 16S ribosomal RNA (partial)
    3     1 15392798 15396046           23S ribosomal RNA
    4     1 15397439 15397917 16S ribosomal RNA (partial)
    5     1 15398097 15398991 16S ribosomal RNA (partial)
    6     1 15399995 15403243           23S ribosomal RNA

I want to plot these features on the karyotype of C. elegans. I am using the chromPlot package from Bioconductor.
I am using the following code:

Subsetting the original data frame to get a data frame for each feature type: 5S, 16S and 23S rRNAs
annot1 <- dplyr::filter(cel, Name == '5S ribosomal RNA')
annot2 <- dplyr::filter(cel, Name == '23S ribosomal RNA')
annot3 <- dplyr::filter(cel, Name == '16S ribosomal RNA (partial)')

and then plot:
chromPlot(gaps = cel, figCols = 2, annot1 = annot1, annot2 = annot2)

This is OK because the two features are far away in the genome. However, if I try:
chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot2,
          annot3 = annot3)

one of the features masks the other because they are close to each other in the genome. How can I make the green bar stack on top of the yellow one instead of covering it?

Did you tried changing the order ? annot1 = annot1, annot2 = annot3, annot3 = annot2) If the annot_n are plotted in that order, simply changing the size will change the results as desired

good idea, but then the green bar is half masked by the yellow bar and this gives the wrong impression of size.

I'm not sure if this is possible at all with chromPlot. Are you sure there is the possibility to make this work as you wish? (have you seen it elsewhere?). Maybe this is a question for the developers themselves

p. 28 of the manual seems to make it work: https://bioconductor.org/packages/devel/bioc/vignettes/chromPlot/inst/doc/chromPlot.pdf

Then you might need to check what "chrSide" parameter does, because is the only difference I see. (And I am not sure that each bar is just added at the top of the other bar, but that could be the effect of chrSide)

chrSide decides on which side this histogram is. I now put 2 and 2 together and solved my problem, thanks.

chromPlot(gaps = cel_gap, figCols = 2, annot1 = annot1, annot2 = annot3, annot3 = annot2,
          chrSide=c(1, 1, -1, 1, -1, 1, -1, 1))

The chrSide argument decides on which side the histogram goes, so it can be used to solve my problem. The bars are not stacked but the result is the same: more accurate visualisation.
Thanks to Llopis for the guidance.


text:
How to concatenate FASTAs with respect to chromosome?
I have the following two small FASTAs, dna1.fa and dna2.fa:
# dna1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG

# dna2.fa
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

If I concatenate these two files using cat, the piped result is as follows:
$ cat dna1.fa dna2.fa > total1.fa
$ cat total1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Here, as expected, these are two files which have concatenated together. It seems this is usually what is done when one concatenates several FASTAs together in bioinformatics.
Question1: Is it possible to concatenate these "based on the chromosome"? Here is the desired result:
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Note that these are exceptionally small FASTAs---usually one could be dealing with files of several GB in size.
Question2: The meta question here is, are there any standard bioinformatic analyses whereby the "concatenated by chromosome" FASTA would be more useful than the straightforward concatenated FASTA via cat alone?

Why would you do this? You're creating contiguous sequences out of (presumably) non-contiguous sequence?

Also you should not assume FASTA format is one sequence per line

What's the actual context here? Are these sequence alignments? Aligned FASTA would not follow the format of your desired output

Using the TblToFasta and FastaToTbl scripts I have posted before, combined with classic *nix utils, you can do:
$ join -a1 -a2 <(FastaToTbl dna1.fa | sort) <(FastaToTbl dna2.fa | sort) |
    sed 's/ /=/2' | TblToFasta | tr = '\n'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation
The FastaToTbl will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to join, which will join its input on lines where the 1st field is the same (by default). The -a1 -a2 means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.
Since the join will add a space between the two joined lines, we deal with that by replacing the second space on each line with a = (sed 's/ /=/2') and then pass the whole thing to TblToFasta to revert to fasta format. Finally, we use tr to replace the = with a newline so you can see where each sequence came from as you show in your example.
For longer sequences where you want to avoid sorting the input, you can do:
$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2
         }
         END{
              for(chr in seq){
                printf ">%s\n%s\n", chr, seq[chr]
              }
         }'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation
Here, the tbl format sequences (see above) are passed to an awk script which does the work:

seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2 : if the 1st field ($1, the sequence identifier) is already in the array seq, then set the value of seq[$1] to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field ($2, the sequence). If it hasn't been seen before, then set the value of seq[$1] to the 2nd field, the sequence.
END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]} : after processing everything, go through each element in the array seq and print out the data in fasta format.

However, that will mean keeping both files in memory which may be a problem for huge sequences. You can avoid that by using temp files:
$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            if(!seen[$1]++){
                print ">"$1 > $1".tmp"
            }
            print $2 >> $1".tmp"
        }' && cat *.tmp > out.fa && rm *.tmp
$ cat out.fa
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation
This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

if(!seen[$1]++){print ">"$1 > $1".tmp"} : if this is the first time we've seen this ID, print it as a fasta header into the file $1".tmp" (e.g. chr1.tmp).
print $2 >> $1".tmp" append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (rm *tmp).

As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.
Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (>chrN) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.
But I don't really know why you would want to do it. You're the one asking for it, after all! :)

Can you answer their question 2, the why?


text:
Relationship between sequencing lane and ngs dataset
I am fairly new to NGS data analysis and I am struggling to understand the exact relationship between a sequencing lane and an NGS dataset. I should add I don't work in the lab, I only do bioinformatics.
I understand the basics of how to make NGS libraries (very briefly: shearing, ligation of adapters, optional PCR, ligation to solid surface, following by sequencing-by-synthesis or other method). What I am struggling with is the visualisation of what happens once you have made your libraries. My questions are:

If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?

If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell? Your provider would then give you raw FASTQ files that have been identified using this barcode? So you would effectively received one (or 2, depending if it's single/paired-end sequencing) FASTQ file per patient right?

If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?
Apologies if these are basic questions.

If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?

Correct, you have to prep each sample into a library.

If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell?

You could, or you could pool it with other samples (assuming they have compatible indices). As @swbarnes2 mentioned, you don't see a lot of samples run on a single lane these days. Sequencing centers like to split things across lanes in case one lane fails for some reason or there are bubbles in the flow cell or other technical issues.

Your provider would then give you raw FASTQ files that have been identified using this barcode?

Typically, yes, and if your sample was run on multiple lanes they'll often be merged for you automatically. This isn't always the case. I've seen a number of people receive samples either still split by lane (annoying, but you can just cat things together) or, much worse, not demultiplexed at all. In the latter case, people often receive 3 (or 4) fastq files per lane, one for each read and another for each index read. It's rare that sequencing facilities do this, but it happens.

If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?

You could sequence a sample per lane. Whether that makes economic sense is up to you to decide.

re but you can just cat things together: do you have any reference for this? I had demultiplexed files, but per lane files were separated. When I done a kmer spectra analysis of individual libraries using KAT, I could find two lanes being far more distant from the other lanes. Is this possible? I would open a new question, but the reference should be in this answer I think.

That's only possible if there are optics or fluidics errors in the sequencer, in which case the run will generally appear lower quality. The sequencing facility would have noticed this if they were paying attention (both the control software and the demultiplexing software give relevant QC metrics). Normally the resulting bases are lower quality then as well.

It is worth noting that the Illumina NovaSeq and NextSeq address all of the lanes of the flow cell from a single loading point (by default). So on those instruments, you end up with all libraries loaded across all lanes. Therefore for each index you will have L001, L002 (+ L003 and L004 if appropriate) FASTQ files (and R1/R2, again if appropriate).
As @Devon mentions, these files can be cat-ted together, but for some processes modelling the lane effect does matter (BQSR, for example).

Indicies matter, lanes do not, unless you have a screw up on the instrument that messes up some lanes but not all.  I'm not sure if anyone skips indices just because they put only one sample on a lane nowadays...it's valuable to be able to use the indices to catch sample mix-ups.


text:
Annotating splice junctions from tophat/STAR output
Is there a way to annotate the splice junctions output from tophat/STAR output?
What I mean by annotate is can I know if it was involved in an alternative splicing event say skipped exon, MXE or retained intron...?
I did some research looks like rMATS does it but we do not have access to the data all we have is output from tophat(.bed file) and STAR(SJ.out) files
Data available:
Splice junction files from star and tophat2 for some reason they deleted the bams in the same directory (from what i know they purged all files with size greater than something)
Species: human aligned to GRCh37(ensembl).
Goal :
Annotate splice junctions based on event it is involved an event skipped exon, MXE or A3ss or A5ss.
I did some more research and saw that ucsc table browser has UCSC Alt events as an option the data looks like this
#bin    chrom   chromStart  chromEnd    name    score   strand
0   chr1    67108547    67109226    atacIntron  0   +
22  chr1    111148925   111149075   altPromoter 0   -
35  chr1    220046661   220292759   strangeSplice   0   +
73  chr1    917423  917573  altPromoter 0   -
75  chr1    2489907 2491261 retainedIntron  0   +
93  chr1    21494532    21496156    strangeSplice   0   -
95  chr1    23720526    23724037    atacIntron  0   -
100 chr1    28835417    28837404    bleedingExon    0   +
153 chr1    84259617    84379035    strangeSplice   0   -

Not sure how they curate this.


text:
Error in seq.default in chromPlot
I am using chromPlot to visualise the genome of C. elegans.
library(chromPlot)
I have created the following data frame with the lengths of C. elegans chromosomes.
Chrom Start      End    Name
1     1     0 15072434 contigs
2     2     0 15279421 contigs
3     3     0 13783801 contigs
4     4     0 17493829 contigs
5     6     0 20924180 contigs
6     7     0 17718942 contigs
7     5     0    13794 contigs

When I use chromPlot(dataframe), I get a plot with 4 chromosomes and the following output/error:
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
Chrom 1 : 15072434 bp
Chrom 2 : 15279421 bp
Chrom 3 : 13783801 bp
Chrom 4 : 17493829 bp
Error in seq.default(minval, maxval, length = howmany) :
  'from' must be of length 1

Does anybody know why I get this error?
I don't get this error when I use the data frame provided by the package with human segments. The structure of the data frame is:
'data.frame':   457 obs. of  4 variables:
 $ Chrom: chr  "1" "1" "1" "1" ...
 $ Start: int  124535434 121535434 3845268 13219912 17125658 29878082 120697156 120936695 121485434 142731022 ...
 $ End  : int  142535434 124535434 3995268 13319912 17175658 30028082 120747156 121086695 121535434 142781022 ...
 $ Name : chr  "heterochromatin" "centromere" "contig" "contig" ...

Compared to my data frame:
'data.frame':   7 obs. of  4 variables:
 $ Chrom: chr  "1" "2" "3" "4" ...
 $ Start: int  0 0 0 0 0 0 0
 $ End  : int  15072434 15279421 13783801 17493829 20924180 17718942 13794
 $ Name : chr  "contigs" "contigs" "contigs" "contigs" ...

I am on MacOS X.
Update
The 20th line of the package-provided hg_gap data frame is:
   Chrom Start   End     Name
20     1     0 10000 telomere

Note that the Start is 0.
I've also tried changing the start from 0 to 1. The error stays the same.
Update 2
Sorting the data frame by chromosome name does not improve the situation.
The traceback of the error is:
6: stop("'from' must be of length 1")
5: seq.default(minval, maxval, length = howmany)
4: seq(minval, maxval, length = howmany)
3: pretty_ticks(minlab, maxlab, 2, ...)
2: draw.scale(y = 0.8 * chr.length[[chrom]] + 0.2 * xylims[3], minval = margin,
       maxval = maxGeneCount, minlab = annot1_plot_range[1], maxlab = annot1_plot_range[2],
       lwd = 4, col = colAnnot1, cex = cex, title = scale.title)
1: chromPlot::chromPlot(gaps)

I have emailed the maintainer of the package Karen Orostica with a link to this post.

I think the error occurs because the start should be 1 or more and not zero. I think the BED format should be 1 based (now it is 0 based).

I've updated my question. I don't think this is the problem based on my experimenting with a 1 start.

Maybe Chrom needs to be sorted, in your example after Chrom 1-4, comes 6 not 5. I am just guessing here.

It seems more like it tries to do something like seq(NULL, to = 50) . If the data is valid, it is a bug on the package. Create an issue as required by the maintainer of the package adding the traceback to see the origin of the error.

Sorting does not make a difference. I have sent an email to the maintainer Karen Orostica.

Are you sure C. elegans has 7 chromosomes? According to ensembl it has only 5 and X.

mitochondrial "chromosome", for my own business I consider the mitochondrial genome like a chromosome but I concede it's technically wrong.

also chromPlot didn't like my character names, but one problem at a time.

I was wondering if that may cause the error. Does it work if you only use the first 4 lines of your df? I mean do you still get the error?

The code works with the first two lines of the data frame but fails if I run it one the first three lines.

This gave me an idea. If I change the chromosome 3 length to 16000000, it works, suggesting it's not allowed to be shorter than the previous chromosome? But ordering by end base did not help.

oooouh, but ordering by end base and removing the mitochondrial genome worked! still not an ideal situation. I wonder if I'm using it wrong or if there is a bug.

@charlesdarwin Maybe that could be written as an answer, just so future reader know how to solve/avoid this problem

@Llopis unfortunately, I could not reproduce it but I've got an alternative that I'll put as an answer.

Using this data frame:
gaps
  Chrom Start      End    Name
1     I     0 15072434 contigs
2    II     0 15279421 contigs
3   III     0 13783801 contigs
4    IV     0 17493829 contigs
5     V     0 20924180 contigs
6     X     0 17718942 contigs
7 MtDNA     0    13794 contigs

and this command:
chromPlot(bands = gaps)

I could obtain this plot:

and then, if you want to add some segments for a subset of bases, you can do:
gaps2 <- gaps[1:2, ]
gaps2$End[2] <- 7000000

chromPlot(bands = gaps, segment = gaps2)

and obtain the following plot:

A more satisfactory answer could put segments of different colour on the chromosome, so bands but this is the best I've got so far.


text:
Meaning of the FORMAT fields of the VCF file coming from GIAB project
After reading the GIAB paper in https://www.biorxiv.org/content/early/2018/05/25/281006 and its Figure 1, I am still having trouble understanding the data inside the GIAB VCF file for HG001 (HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_all.vcf). In particular, I need to understand the information that led to a FILTER value of PASS for some variants.
The first case that I would like to discuss is the following variant:
chr1    110176467   .   A   G   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable    GT:PS:DP:ADALL:AD:GQ    0/1:.:512:114,106:145,148:760

Translating the information in the line to English, this is what I get:

We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
However, the callsets from CS_HiSeqPE300xGATK,CS_CGnormal, and CS_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
Across all the platforms, DP=512
ADALL (Values 114,106) is coming from “all" the datasets. Are those “all” the set (HiSeqPE300x + Cgnormal + 10XChromium), or the set (Illumina + CG + 10X + IonExome + SolidPE50x50bp + SolidSE75bp).
AD (Values 145,148) is coming from the same “all” datasets as ADALL, but they are unfiltered, so 145>114 and 148>106.

All this make sense, but I would like to know what case of Figure 1c was applied to arrive at the PASS value.
The second variant that I would like to mention is this one:
chr1    5705293 .   T   C   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_HiSeqPE300xGATK_filt   GT:PS:DP:ADALL:AD:GQ    0/1:.:599:121,139:49,28:627

Again, my translation would be:

We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
However, the callsets from CS_HiSeqPE300xGATK,CS_CGnormal, and CS_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
Also, the data from CS_HiSeqPE300xGATK_filt was filtered/discarded
Across all the platforms, DP=599
ADALL (Values 121,139) is coming from “all" the datasets. I assume (Cgnormal + 10XChromium), because CS_HiSeqPE300xGATK_filt was filtered.
AD (Values 49,28) are unfiltered.

But notice that 49<121 and 28<139. How come that there are less unfiltered reads (AD)? Should not be the case that AD > ADALL for all instances?
I don’t know how to reconcile the two variants. Aren’t they providing opposite information in the AD and ADALL fields? Which specific datasets and callsets are involved in getting AD and ADALL for each variant?
I would also appreciate answers to a couple of other questions:

Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?
Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

For reference, here are the descriptions in the header of the VCF file that I have been using:
##INFO=<ID=platforms,Number=1,Type=Integer,Description="Number of different platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=platformnames,Number=.,Type=String,Description="Names of platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasets,Number=1,Type=Integer,Description="Number of different datasets for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasetnames,Number=.,Type=String,Description="Names of datasets for which at least one callset called this genotype, whether filtered or not">
##FORMAT=<ID=ADALL,Number=R,Type=Integer,Description="Net allele depths across all datasets">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Net allele depths across all unfiltered datasets with called genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Total read depth summed across all datasets, excluding MQ0 reads">
##INFO=<ID=datasetsmissingcall,Number=.,Type=Integer,Description="Names of datasets that are missing a call or have an incorrect call at this location, and the high-confidence call is a variant">
##INFO=<ID=callable,Number=.,Type=String,Description="List of callsets that had this call in a region with low coverage of high MQ reads.">
##INFO=<ID=filt,Number=.,Type=String,Description="List of callsets that had this call filtered.">
##INFO=<ID=arbitrated,Number=1,Type=String,Description="TRUE if callsets had discordant calls so that arbitration was needed.">

The filtering criteria should be defined in the header of the VCF file. Can you include the header from your VCF file?
VCF Format and example header:
FILTERs that have been applied to the data should be described as follows:

##FILTER=<ID=ID,Description=”description”>

header:
##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">

Thanks for the help. I posted the header fields that I am having trouble with.

For your second question, I don't think that variant matches any of the examples shown in Figure 1c. Note that that panel's caption says it shows four examples of arbitration. Per you last sub-question, I don't think this variant was subjected to arbitration.
For your second question, I believe the subtlety is in the field descriptions:
ADALL = Net allele depths across all datasets
AD = Net allele depths across all unfiltered datasets
with called genotype
My understanding is that ADALL is a count of allele depths that includes samples where the variant was filtered out or the genotype was not concordant; something to that effect. So generally I would expect ADALL > AD.
I'm honestly not sure exactly what they mean by "unfiltered datasets". That could mean "counting all reads without filtering on quality" or "callsets to which variant filtering was not yet applied". Maybe there's another subtlety there that explains why the two variants have opposite AD vs ADALL proportions.
For your additional questions:

Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?

According to the paper, having the variant in low-coverage regions isn't disqualifying if the genotype calls are concordant, which makes sense to me. Low coverage makes me more suspicious, but if multiple callsets are in agreement, I'm still willing to consider the call good (unless of course they all exhibit the same error mode, which is the main potential weakness of the consensus approach used by GiaB).

Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

Right, if that flag is not present, it means arbitration was not needed/applied.

Thanks. I also think that there is a subtlety somewhere. Logically it does not make sense that ADALL >  AD, no matter if called or not.


text:
How to download dbSNP database?
I am bench-marking the performance of ADAM/Avocado on NGS whole-genome sequencing datasets. I am trying to repeat the tests reported in Frank Nofthat's (author of ADAM and Avocado) thesis, however, I am running into an issue acquiring the dbSNP database to use for Base Quality Score Recalibration.
The Base Quality Score Recalibration command (page 128-129):
adam-submit -- transformAlignments aln.adam bqsr.adam -recalibrate_base_qualities -known_snps dbsnp.adam

requires dbsnp.adam, which is described as "the dnSNP database" (converted to ADAM format) with reference to the NCBI database of genetic variation. This "database" can be found at this FTP site however I do not know which parts of these files to download and convert to ADAM in order to perform this step of the transformation pipeline.
Also, in order to use this genetic variation databaes do I need to make sure that the reference genome of "known variants" is the same genome that the sequencing data was aligned to?
Thank you.

In the dbSNP FTP site you linked, you need to go into organisms and select your organism of interest (human obviously). You can then select which release of the SNP database you want (eg. 150 vs 151) [the modified dates hint at when they were released] and which genome build you want (eg. hg37 vs hg38).
Below I am looking at the .vcf files. I assume this what you want because they can be converted to ADAM format via the adam vcf2adam command.
For example human (newest genome and newest dbSNP release):
ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/
Your next decision is whether or not you want to use ALL the annotated SNPs in the database or just the common ones. For the former you choose the All_20180418 file and for the later the common_all_20180418 file. I'm pretty sure the convention here is that the date in the filename is when it was uploaded/updated.
In your case I'm not sure what makes the most sense. Firstly, are you trying to replicate the initial conditions of when that software was published or do you want simply to benchmark the method with the best available data. Besides that, I imagine using all them is computationally much more taxing and if you are just coming up with general trend then using the common ones seems viable.

This is super helpful. One more question though: does my input BAM file need to have been aligned to the same reference genome (in this case GRCh38)? If so I will need to re-align my input BAM to the CRCh38p reference genome before performing this step.

Yes ideally although it looks like they also have an hg37/hg19 version


text:
How to convert files to ADAM format?
I would like to convert BAM and VCF files to ADAM format.
How do I do that?

These transformations are quite simple with the adam-submit script packaged with ADAM.
Transformation of BAM files
adam-submit -- transformAlignments sample.bam sample.alignments.adam

Transformation of VCF files
adam-submit -- transformVariants sample.vcf.gz sample.variants.adam

Could you provide a link to the package? It seems like it is this repository but I'm not sure.


text:
How do I create a for loop to filter through different FDR values?
This is probably a really quick/simple fix but I am a Noob to R programming/programming in general. I am trying to create a for loop that runs the test multiple times with different FDR values (5,10,20). How do I make it run 3 times with different values and then produce 3 different datasets with names like... atac.glmtop.males.opening.fdr5, atac.glmtop.males.opening.fdr10, atac.glmtop.males.opening.fdr20? Also is there a way to make the write.table function work so that it exports files with the fdr5, fdr10, fdr20 in the name?
fdr.nums = c("5","10","20")
fdr.nums.vector = c(5,10,20)

for (i in 1:length(fdr.list)){
  atac.glmtop.males.opening.fdr[i] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr.vector[i]/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[i]

  write.table(
    atac.glmtop.males.fdr[i].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",i, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[i].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ i + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

I'm voting to close this question as off-topic because it is not about bioinformatics.

For loops in most languages can be done in a variety of ways. What you have currently coded, loops over values 1 to 3 by using
for (i in 1:length(fdr.list)){
...
}

However, you probably want to use a different form and loop over each value in that vector using something more like:
for (fdr in fdr.nums.vector){
...
}

This way you have convenient access to the FDR value in each loop via the fdr variable. Which you can then use in constructing your output filenames. Thus your code would become something like:
fdr.nums.vector = c(5,10,20)

for (fdr in fdr.nums.vector){
  atac.glmtop.males.opening.fdr[fdr] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[fdr]

  write.table(
    atac.glmtop.males.fdr[fdr].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",fdr, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[fdr].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ fdr + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

Advancing From for Loops to apply
As a beginner, for loops seem quite logical and easy to understand. However, R has a more powerful way to achieve the same thing through what is called "vectorisation" using the apply function or one of it's relatives lapply, sapply etc. Once you have understood and mastered for loops, you should take a look at apply, especially if you ever need to do a for loop many times or you find yourself doing nested for loops. R is notoriously slow with for loops and apply will likely be many times faster (potentially orders of magnitude). Here are some links for when you are ready to advance from for loops to vectorisation.

https://www.r-bloggers.com/how-to-write-the-first-for-loop-in-r/
https://www.r-bloggers.com/using-apply-sapply-lapply-in-r/

You are right about vectorisation, however, apply is not vectorisation but uses for loops, see this for details: https://www.r-bloggers.com/vectorisation-is-your-best-friend-replacing-many-elements-in-a-character-vector/
or the document 'The R Inferno'.

“As a beginner, for loops seem quite logical and easy to understand” — this is really only true for people who are used to imperative programming. Properly taught, a functional-first approach is simpler for beginners. I strongly recommend skipping for loops in teaching and learning (or at least postponing them until much later in a course). Especially since, in R, they require all kinds of special explanations.


text:
Find the date that sequencing was performed
Background I've obtained RNAseq datasets from ENCODE and I'm looking to identify external variables contributing to sequence variance among these samples.
External variables For example, clustering effects on the date the sequencing was performed.
Questions How do I find out the precise date which each sample went through the sequencing machine?

ENCODE databases are here on NCBI
The ENCODE website is here

What kind of data do you have? Do you have the count table or the original fastq files? Where did you get this files from (If you got it from someone could you ask him/her) ?

I downloaded them directly from the ENCODE website. I thought to ask here before e-mailing.

And why do you need the date? It is for batch finding? I don't know the metadata for the project but usually there documentation about how is this encoded

@Llopis: Yes, I'm looking to correct for variance associated with dates. Maybe I'm stretching it, but I'm really interested in parsing through the noise of the different datasets at all factor levels

Will you want to know which person did each part of the process? Or in which center, machine and lane it was done? I agree that this should be an objective, but I think it is unreasonable to ask this of this big projects, any way you can ask them. But I see there is a "Date released" filter value, so you can use that as a proxy

Yes, I agree-- for some parts it may be overkill. Date released corresponds to when the consortium released the data to the public. In fact, the sequencing lane and machines are in fact available in the library metadata! I thought maybe there is a way to get date information from a fastq file... I've read that dates are fingerprinted onto microarrays.


text:
Genotyping with Illumina HumanOmni1-Quad before whole genome sequencing
I read a paper about  (whole genome sequencing analysis of 100 Southeast Asian Malays. It says "Prior to library preparation, each sample was genotyped on the Illumina HumanOmni1-Quad as an initial DNA QC and for in-house assessment of gender and familial relations". I can hardly understand how Illumina HumanOmni1-Quad worked here for initial DNA QC and for in-house assessment of gender and familial relations. Could you please kindly explain to me a bit about this sentence?
Many thank in advance

The Illumina HumanOmni1-Quad beadchip is a microarray device consisting of 1,140,419 markers which have derived from the 1,000 genomes project. The markers chosen are is high-value regions of the genome: ABO blood typing SNPs, cSNPs, disease-associated SNPs, eSNPs, SNPs in mRNA splice sites, ADME genes, AIMs, HLA complexes, indels, introns, MHC regions, miRNA binding sites, mitochondrial DNA, PAR, promoter regions, and Y-chromosome.
The markers come together to form a sort of individual fingerprint. Germline SNPs are heritable and so can be traced up the line. Some SNPs are more prominent in subpopulations, something that Ancestry.com and other similar services leverage to predict ethnicity.
This paper utilized genomic DNA extracted from the Singapore BioBank. BioBanks tend to record age, lifestyle, family history, and health information alongside the genomic DNA. I (conservatively) assume that the Singapore BioBank at least has familial relations and gender information tied to the DNA. My interpretation is that the ivestigators SNP genotyped using the Illumina Beadchip the genomic DNA from the bank and compared their results with the "ground truth" provided by the genetic bank.
Illumina brochure for the HumanOmni1-Quad beadchip.


text:
High Phred Quality score VCF file
I am working on ION Torrent S5 XL and I am facing to an issue with VCF output generated by the variant caller of thermofisher. It stating on VCF 4.1 that Phred Quality Score is on a range of 0 to 255. In my VCF I have values over than 255. Do you know if there is a way to transform these values?

Let me rephrase to see if I understood: your output exceeds the max value explained by its manual? Could you post an example? My only guess is that either the structure is wrong or there is a bug, or the program is not well documented.

Hi Llopis, here's an example: chr17 7577104 . AGGA AAAA,AAGA,AGA,AGAA,AGGG,AGTA 35144.8 PASS. 35144.8 is QUAL Score into the VCF

Please can you [edit] the question to include it in the right formatting? Also could you post the command you used to get this VCF ?

See this post on the GATK forums.
It is normal to see values above 255.  The QUAL score is the posterior probability that all samples in your callset are homozygous reference. I quote Sheila's response from the thread posted above:

The QUAL score is the Phred-scaled posterior of AC = 0. We use the AC priors and the PLs to get the likelihood of the data given each AC, then use those to get the posterior probability for each AC. From there, the calculation is 1 - Pr{AC > 0}.

Please see the GATK post in Methods and Algorithms for the nitty gritty details. Beware, Bayesian statistics abound.

GATK's HaplotypeCaller and freebayes are two more variant caller that output QUAL values greater than 255. The reason for that is, that each read have its part on calculation that value. This results in the higher the read depth at this position the higher the QUAL value.
To get something useful out of this value, the usual thing that is done, is to calculate the quality per depth (QD) by just divide the QUAL value by the read depth.
Note whenever you try a different variant caller, you have to find out whether the QUAL value is read depth dependend or not.


text:
Retrieving NCBI Taxa IDs from refseq or GenBank assembly accession
I have about 10,000 genome files all named by either refseq or genbank accession number, do you know if it's possible to convert these numbers to the corresponding NCBI taxon ID or species?
for example:
GCA_000005845.2 to 79781

In the case of E.coli
Edit:
The file names look like this:
GCF_000337855.1_ASM33785v1_protein.faa.gz
GCF_000091125.1_ASM9112v1_protein.faa.gz
GCF_000184535.1_ASM18453v1_protein.faa.gz

My operating system is Ubuntu 16.04

Hi and welcome to the site. Please [edit] your question and show us a few examples of the actual file names. If you have both refseq and genbank accessions, show examples of both. Also, please mention your operating system since some solutions might depend on it.

Thanks for the edit. So, are your filenames always accession.versionassemblyprotein.faa.gz?

Yes, although they are all in different folders, for example: /refseq/archaea/GCF_000337855.1/GCF_000337855.1_ASM33785v1_protein.faa.gz

Turns out I'd already written some code that did it, human memory is a funny thing.
This uses biopython to split the field description to where the species is. May not work for all NCBI files, but seems to work on most.
import Bio
from Bio import SeqIO
from Bio import AlignIO

for record in SeqIO.parse (FILE, "fasta"):
    Speciesname = record.description.split('[', 1)[1].split(']', 1)[0]

You can use eutils to get this information as shown below:
esearch -db assembly -q 'GCF_000337855.1' | esummary | xtract -pattern DocumentSummary -element AssemblyAccession,Taxid
GCF_000337855.1 469382

As the question mentions python, here is a python example using biopython.
import json
from Bio import Entrez

def accession2taxid(acc: str, db="nucleotide") -> str:
    handle = Entrez.esearch(db=db, term=acc)
    record = Entrez.read(handle)
    gi = record["IdList"][0]
    handle = Entrez.esummary(db=db, id=gi, retmode="json")
    result = json.load(handle)["result"]
    taxid = result[gi]["taxid"]
    return str(taxid)

So if you want to use an assembly accession, as in your example
>>> accession2taxid("GCA_000005845.2", db="assembly")
'511145'

In your example, you said that GCA_000005845.2 should return 79781, but that Taxon ID is for "uncultured archaeon AOC672". The taxid in the example I just showed returns "Escherichia coli str. K-12 substr. MG1655", which is what you want.
Your question also mentions GenBank accessions, so here is another example
>>> accession2taxid("CU458896.1", db="nucleotide")
'561007'


text:
Reduce number of transcripts in a highly variable de novo transcriptome assembly
I have a de novo assembly using both multiple SRA and locally sequenced transcriptomes. I started with 270M PE reads from 9 tissues. Here are the assembly stats generated with TrinityStats.pl:
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  543668
Total trinity transcripts:      1111089
Percent GC: 33.77

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2117
        Contig N20: 1324
        Contig N30: 945
        Contig N40: 714
        Contig N50: 561

        Median contig length: 344
        Average contig: 496.17
        Total assembled bases: 551293315

#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2172
        Contig N20: 1333
        Contig N30: 925
        Contig N40: 681
        Contig N50: 523

        Median contig length: 323
        Average contig: 473.63
        Total assembled bases: 257499729

As you can see, the total Trinity 'genes' number is very high. I suspect that this is due to a high biological variability of the transcriptomes (coming from organisms sampled in Italy, China, Spain...).

How can I assess whether the problem is really the biological variability?
How can I reduce this number? I'd rather not downsample the reads. It looks like CD-Hit might do what I'm looking for...

As you've suggested, CD-HIT works for reducing transcript numbers. We used a mixture of expression-based filtering and CD-HIT for reducing transcript counts for our genome-guided transcriptome assembly. This reduced numbers by a lot, without much change in BUSCO scores:

Map RNA-seq reads to Trinity-generated transcripts using Salmon
Use the expression of BUSCO genes to set a credible signal cutoff
Subset transcripts based on this threshold [for us it was 50 counts]
Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
Run cdhit to collapse similar transcripts

Here's our CD-HIT command line:
cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta

More details here:
https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9

Thank you for your answer. Prior to using CD-HIT, I decided to filter the reads according to BLAST hits with a custom pipeline. Only 150k transcripts had a hit, and of these, only 56K had a hit with the suitable taxon. This means that almost a million transcripts don't match with anything, and probably CD-HIT won't help much. Any idea on why this is happening? Also, given that BUSCO retrieved 975/978 core genes (so my assembly is basically complete), should I care about this number?

CD-HIT does an internal match; it doesn't care about matches to external databases. That BUSCO score is great; even C. elegans isn't 100%.

Yes, but my question is: should I worry about the fact that I have >1M transcripts, and of these only 5% seem to correspond to a mRNA? What's in the other 950k transcripts reconstructed by trinity?

Ah, and just for the record, today I ran CDHIT and it obtained 780k clusters, so yes, it improved the situation, but not to a great extent. This still looks unrealistic to me.

To better answer that question, I need more information about what you were sequencing, how the sequencing was done, and what sequences are coming up. For example, if there's a high chance of bacterial contamination and sequencing is for total RNA including microRNAs, then 780k is something that could happen.

Mind if we talk about it in chat? I'm afraid this is getting long for comments :)

Seems to me that your "genes" are very short. What sort organism (prokaryote, lower eukaryote, higher eukaryote) is this? I would look to see if your "genes" contained ORFs, and only keep those with sufficicntly long CDS. You will miss uORFs and ORF fragments this way, but you'd have much more confidence in what you did find.

It's a clam. Aren't these values OK for a transcriptome? Where can I find good reference values?

Higher eukaryote? Talking about giraffes? http://scienceblogs.com/evolgen/2006/10/18/do-higher-eukaryotes-smoke-mor/


text:
Negative scale-free topology
I am trying to build gene co-expression networks with WGCNA by combining in-house and publicly available RNA-seq datasets. I am interested in identifying gene networks associated with different environmental conditions. I used SVA to adjust for batch effects and "protected" different environmental conditions in the model as I don't want to lose information about them. However, when I used the "batch-adjusted" data to create networks, all the scale-free topology indices were negative (max = -0.05).
I tried generating the networks with softpower = 16 as it is suggested for sample sizes 20-30. However, only five networks were generated which is far less than I expected.
What might be the reason for getting negative scale-free toplogy indices?
Can combining studies from different environmental conditions cause the networks to be not biologically meaningful?

Welcome to the site. I have some questions before being able to provide an answer: Could you show (the code) how did you modify your data to remove the batch effects with sva? How many samples and genes do you have? Are all the RNA-seqs built against the same genome version? How did you call the network construction, full automatic or manual?

Whenever you have batch or other nuisance variables driving variation you're going to get odd networks. You tried using SVA, which is a good method to use, but it seems in this case that it's able to optimally handle this dataset. I also strongly suggest that you use pickSoftThreshold() to choose a softpower, since that may produce better results (though since you already have negative values I'm not holding by breath). In general, you're going to need to play with either SVA or filtering to remove the batch-related genes that are causing problems.

Thank you for the response, Devon Ryan. Removing the low expressed genes (log2(FPKM) < 1) before using SVA has partly resolved the issue. I still don't get a high scale-free topology, but now I have positive numbers (max = 0.6). The resulting modules look trustworthy.

@lence any update on the question? (you seem the same user, but your user id is different, do you need help to merge both accounts ?)


text:
Calculating L2 norm (Euclidean norm) to set of genes in time series
I have gene expression RNA-seq data of 9 time points (0 hour, 2 hour, ..., 16 hour). I need a plot on which x axis is time point and each data point at a certain time point represents the L2 norm (Euclidean norm) of the difference between the gene expression vector at that time point and the previous time point.
This is the link of my expression data
https://drive.google.com/file/d/1crJSqVdGujy1tFyeoV6E7e-xO9SXteJZ/view?usp=sharing

If you read your matrix into R and saved it as m:
diffMatrix = m[,2:9] - m[,1:8]  # The change per-timepoint
l2n = apply(diffMatrix, 2, function(x) sqrt(sum(x^2)))  # The L2-norm

R also has a norm() function, which can directly produce the L2-norm, but it turns out that this is slower than the simple method shown here. If you click through to the other question you can find somewhat faster implementations of this, but anything will suffice for small datasets like yours.

Sorry, How can I normalise these 8 points to zero time point?     12n= Named num [1:8] 313698 168553 360237 487963

That depends entirely on how you want the normalization done.

Sorry I mean how I can normalize these 9 differences you helped me to calulate  to expression at 0 h, which corresponds to cells growing in ground state conditions?? Because they are very big and I can not plot them

You wouldn't, it doesn't make sense to normalize a difference between two time points to a third one. You could normalize everything before subtraction though.

Right, actually the goal was showing global gene expression dynamics as measured by bulk RNA-seq. I was going to plot a data point at a certain time point represents the 2-norm of the difference between the gene expression vector at that time point and the previous time point. However these differences are too big to be plotted. Thank you

Try something like plotting log2(TPM) in violin or joy plots. If there are changes in the distribution of expressed genes over time then that should make it apparent.


text:
No counts for added gene in cellranger (scRNA-seq)
I have a set of scRNA-seq samples enriched with FACS for cells expressing a specific gene reporter (TdTomato). In particular the gene I want to report has positive counts in the resulting matrix for 97% of cells.
I followed CellRanger documentation to create a custom reference. I added TdTomato on an additional chromosome and I added an entry in the gtf:
TdTomato    .   exon    1   1431    .   -   .   gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato"

The command I'm using is:
cellranger-2.1.1/cellranger mkref --genome=refdata-cellranger-mm10-2.1.0.tdtomato --fasta=genome.tdtomato.fa --genes=genes.tdtomato.gtf

where genome.tdtomato.fa is the original reference genome (mm10, as provided by cellranger) with an additional chromosome at the end with header:
>TdTomato

and genes.tdtomato.gtf the original reference genome gtf with the additional line indicated previously at the end of the file.
I actually tried both on the '+' and on the '-' strand, however I guess this should not affect the read mapping, since TdTomato is on an independent chromosome.
I ran cellranger count using the custom reference:
s="SAMPLE"
cellranger-2.1.1/cellranger count --id="$s"_tt --transcriptome=refdata-cellranger-mm10-2.1.0.tdtomato --fastqs=original_fastqs_directory/"$s"/fastqs --sample="$s"

I am using the default parameters of CellRanger.
However, I found that less than 1% of cells seems to express TdTomato. If I use the '+' strand I get slightly more counts (2.4% of cells have positive TdTomato).
What I checked:

With IGV I can see that there are reads mapping to TdTomato
Using samtools I am able to see the same number of reads mapped on TdTom either with the + or - strand
However, following this documentation I noticed that by executing the CellRanger filtering steps, more reads are filtered in the case of the - strand

This reduces the number of reads used for counts to 44 (2.4% of the original 1807 reads, - strand), 240 (13.3%, + strand) for TdTomato, 71620 (21.2%) for the gene reported by TdTomato. Interestingly, there are fewer reads for TdTom than counts, so there may be some kind of normalization afterwards?
The gene reported by TdTomato is longer, 11394 bp than TdTom, 1431 bp. The total length of the exons for the two genes is 5349 vs TdTom: 1431. Does this influence exponentially the amount of reads mapped to the two genes?

Did you check the bam file for reads mapping to TdTomato? If you believe it is a mapping issue you can always run HTseq count on the cell ranger bamfile. It would help if you add the actual input command for cellranger mkref and cellranger count to trouble shoot your problem

As an alternative you can make a reference for an alternative mapper such as  bwa or STAR, count the mapped reads and see how this corresponds to the reads found in the cellranger bamfile. If there is a good correspondence you most likely did not enrich, or used the wrong sequence for your reference.

The strand will not affect read mapping, but will affect read counting. How many TdT counts do you get with different strands? If the counts are similar, you are probably not capturing the real TdT.

@Mack123456 I've added the issued input command, thanks. I'll also look at the bam file with HTseq.

@burger The count is similar, more or less 1% of cells with >0 counts and maximum of 2 or 3 counts.

10X RNA-seq tends to only reliably capture the most abundant transcripts expressed in a cell.  The TdT gene may not be expressed at a level high enough relative to the other transcripts to be reliably captured in your experiment. You may be able to identify additional transcripts if you increase the sequencing saturation.

@gc5 Do you have any update regarding this?  I have done two experiments on cells sorted by expression of tdTomato, one with 10X and another using DropSeq, with a custom tdTomato reference and in both experiments the number of tdTomato transcripts is very small (1/2000 cells in the 10X experiment, and 91/900 in the DropSeq experiment).  I can't figure out if it's a problem with my reference or a biological/technical problem on the sequencing side.

@wflynny unfortunately not yet. I'm still checking I haven't done some mistakes during mapping, but the steps are pretty much consolidated so I don't think there's much room for mistakes. My guess is that TdTomato reads are filtered by cellranger during count, but I don't know why yet. Please update here if you find a possible solution. Thanks

@gc5 With samtools idxstats possorted_genome_bam.bam | cut -f 1,3 for my 10X experiment, I get 130 reads mapping to my "tdTomato chromosome" so I know there are at least some reads mapping.  Very puzzling.

@wflynny I've updated my question with more details, hope this is interesting for your case

@GWW I am thinking that your comment is on point, maybe there are very few reads for TdTomato to begin with and 10X scRNA-seq cannot reliably capture them. Is it basically dropout or another type of issue? Can you provide me some articles to read about this problem? Thanks

The first column of your GTF must correspond exactly to the name of the sequence in your FASTA file. Typically this is the chromosome so a GTF line:
chr1    hg38_knownCanonical exon    11106535    11262507    0.000000    .   .   gene_id "gene1"; transcript_id "tx1";

This must have a corresponding sequence in the FASTA file:
>chr1
ATCG ATCG ATCG ATCG ATCG ATCG...

Therefore you need to add an additional sequence to your FASTA file with >TdTomato. This must have exactly the same sequence as the construct that you are expecting reads for. The run cellranger mkref and cellranger count on these FASTA and GTF files. Cellranger also can take a GTF with “CDS” as a feature in addition to “exons”. The STAR alignment algorithm used by cellranger is splice aware so it can handle split reads but you will need to ensure the GTF corresponds to exons of the mRNA/cDNA rather than the full sequence.
10x Genomics is a platform with a high dropout rate, especially for lowly expressed genes. You should bear this in mind for what you expect to observe. It may be worth comparing the mapping rate of your construct to housekeeper genes to see if this is the reason for the low counts. Since whole mRNA molecules are captured on the beads before preparing the library, there should not be a difference in UMI counts between long and short genes.


text:
Why is FPKM still used for gene expression studies?
Perhaps my understanding is wildly misguided, but I'm seriously confused about why people are still interested using FPKM values for cross sample gene expression analysis. My understanding is, that for some gene ${ i }$, its $FPKM$ is calculated as follows:
$$ FPKM_i \propto \frac{ \mathrm{fragments\ mapping \ to\ gene \ }i } { \sum_j \mathrm{ fragments\ mapping\ to \ gene\  } j } \bigg/
\mathrm{effective\_length}(i) $$
For comparing gene expression within a single sample, this doesn't seem too harmful to me. In fact, this seems like a good idea.
However, the use of $FPKM$ or even the famous $log_2(FPKM+1)$ in a similarity analysis seems to me plain wrong. True differences in RNA abundances of a couple of genes would significantly impact the FPKM values of the other genes leaving  perhaps confused about the true differences between the samples.
Am I missing something here? Is it ever meaningful to use FPKM values for analyzing across samples?

Is it ever meaningful to use FPKM values for analyzing across samples?

One should never use FPKMs for anything important. They can occasionally be useful for plotting, but even in that case one needs to construct the FPKMs from properly normalized data and not use the original definition.
It's becoming ever rarer to see people using FPKMs these days, thankfully.

Thanks, I also saw your response here https://www.biostars.org/p/317417/ as the answer to my followup.

So what is a proper alternative? TMM?

@LinuxBlanket TMM or another robust normalization method, yes.


text:
Which secondary structure prediction algorithms don't query the PDB/DSSP for known homologous sequences?
There are several secondary structure prediction algorithms. I am interested in the algorithms which predict secondary structure without accessing the PDB to cheat by using the DSSP or STRIDE entry for known homologous structures.
Since DSSP and STRIDE are not performing predictions, rather they analyse the PDB structure to determine whether each residue is a helix, strand or coil.  This in turn means that, the 'so called' secondary structure prediction algorithm, also, is not performing predictions (except, of course, in the cases where no homologous structure/sequence is found in the PDB).
Which secondary structure prediction algorithms don't cheat by using DSSP or STRIDE for homologous structures?
Here is a list of secondary structure prediction algorithms, which are accessible via a consensus secondary structure prediction server called @NPS (currently down: https://npsa-prabi.ibcp.fr/NPSA/npsa_seccons.html)
DPM
DSC
GOR1
GOR3
HNN
MPSA
PHD
PREDA
SOPM

Have you tried googling for ab initio secondary structure prediction?

@DevonRyan It is a good idea, and doing so does provide some useful information.  However, ab initio modelling is not what I am interested in.  Use of knowledge is fine, only it should be generalised, not a specific query to the PDB and DSSP for a specific sequence - in such case, no prediction is performed, only cheating i.e. an ~100% accurate answer.

If using knowledge is fine the querying PDB and such when possible is in no way cheating.

Yes it is cheating.  Not performing a prediction = cheating.

You have a warped definition of cheating if it includes using the most reliable information. Since you implicitly don't want the most likely results, what is your end goal with all of this?

A "predictor" not making a "prediction" is cheating.  I am conducting some research into secondary structure predictions comparing their predictive results.


text:
Any input as to what to study in order to work into genetic engineering from math background?
Any input as to what to study in order to work into genetic engineering from math background?
I have troubles understanding what the "interface" would be between mathematicians and biologists/biotechnologists.
Wikipedia links the first model of an organism having built in 2011:
https://en.wikipedia.org/wiki/Computational_genomics#First_Computer_Model_of_an_Organism
And it based on statistically analyzing biology/chemistry papers. This reminds of a typical "stats research". So there's not necessarily direct interaction between the empirical and the formal. However, for new applications, it's given that there must be some interaction, in order to synthesize knowledge?

I don't understand the question. Do you seek advice to how to proceed in your career or do you want to know how does an "empirical" and a a "formal" science interact?


text:
Batch detection of CRISP proteins in fasta file
Probably a naive question. I am inexperienced.
I am interested in identifying potential CRISP (Cysteine-rich secretory proteins) in a certain tissue transcriptome (ca. 20k sequences in fasta). I have detected signalP and estimated % of cysteine in sequences. However CRISP rely on a certain pattern, and I am not aware of batch-detected algorithms.
Please would anyone know the best/simplest way of detecting CRISP candidates?
[EDIT]
I have been asked for the pattern. I summarise what the literature says:

Several of them have a CAP domain;
There should be a region called 'cysteine-rich domain' (CRD) upon the carboxyl terminal half of the protein containing 10 of roughly 16 conserved cysteines.

I am particularly interested in venom-derived CRISPs, but I am not sure whether they cluster in a certain CRISP family.
References:
Roberts et al. (2007) Structure and function of epididymal protein cysteine‐rich secretory protein‐1.
Yamazaki & Morita (2004) Structure and function of snake venom cysteine-rich secretory proteins.

Could you [edit] your question and tell us what this pattern is? Is it a specific sequence pattern? If you can define the pattern, we could be able to figure something out even if there aren't any existing dedicated tools for this.

@terdon I have read some references and hopefully the pattern described in my edit update comprises what I am looking for. Thanks

Is your goal to detect the actual proteins from the sequences? If you'll settle for detection of transcripts from the genes encoding these proteins as opposed to detection of the proteins themselves, here is an approach.
Aligning your transcriptome to a reference genome (probably with STAR) will be the first step. You have lots of short reads in your FASTA, but in their current state they are no more informative than a box of puzzle pieces. Aligning your reads to the reference will allow you to extrapolate which genes that you short reads correspond with. From here you can quantify the features (genes) in your transcriptome.
Lucky for you, there are only 3 genes in the cysteine rich secretory protein family, CRISP1, CRISP2, and CRISP3. Alignment with a reference (and annotations) and quantification of genes that your reads align to will give you all you need. You'll just have to find the three genes in your final counts table.
I recommend STAR as an aligner with the --quantMode GeneCounts argument. It's very quick if you've got the resources and is geared towards RNA-seq reads.

My fasta is the result of de novo assembly, plus we have a reference genome from a proximate species for comparison. Problem is, these are not intensively-annotated organisms (e.g. not Drosophila nor Mus musculus) thus I am trying to fish out as many annotations of novel transcripts as possible.


text:
Correcting for noise in RT-qPCR gene expression data
I have a training set of RT-qPCR gene expression data (not run in triplicate) for a batch of samples with two phenotypes $A$ and $B$ on which I've trained a "logistic regression classifier".
I also have another smaller set of samples which have been run twice as technical replicates on the RT-qPCR machine in order to study the effect of qPCR noise on my classifier.
I have a particular cut-off probability $P$ for my logistic regression model, such that we must be at least $P$ confident to classify a test sample as $A$ (similarly for $B$), otherwise no classification is returned.
I would like to incorporate the effect of RT-qPCR noise into my logistic classifier in order to return a more accurate classification probability for a particular test sample.
My model for the RT-qPCR noise for a particular sample $i$ is given by
$$X_i\sim N(\mu_i, \sigma^2I).$$
The logit is then given by
$$t_i = \beta_0 + \beta^TX_i,$$
with distribution
$$t_i\sim N(\beta_0+\beta^T\mu_i, \sigma_t^2),$$
where $\sigma_t^2=\|\beta\|_2^2\sigma^2$.
Since $\mu_i$ is a nuisance parameter, we eliminate it by taking the difference of the two replicates for a particular sample
$$Z_i=t_i^{(1)} - t_i^{(2)} \sim N(0, 2\sigma_t^2).$$
We can then form the MLE of $\sigma_t^2$ by computing
$$\hat{\sigma}_t^2=\frac{1}{2N}\sum_{i=1}^NZ_i^2.$$
Therefore, for a particular test sample $x_i$, we compute the expected value of the logit-normal distribution with location $t_i$ and squared scale $\hat{\sigma}_t^2$
$$p_i=\int_0^1x\frac{1}{\sqrt{2\pi\hat{\sigma}_t^2}}\frac{1}{x(1-x)}e^{-\frac{(\text{logit}(x)-t_i)^2}{2\hat{\sigma}_t^2}}dx.$$
The resulting $p_i$ is then our corrected probability of phenotype $A$.
Is this a sensible way to achieve my stated goal?
Here is an example plot of the necessary correction when the squared scale is $\sigma^2=2$,

The black line corresponds to a qPCR assay with no noise.

The violet line indicates the corrected probability.

For example, if the regression model says that a sample has a 75% chance of being phenotype $A$, then in the corrected model, this actually corresponds to only having about a 65% chance of being phenotype $A$.

You're actually not correct about the functional form of the logit. It's nonlinear, but you've written a linear function.
That matters here because once you understand the functional form of logistic regression, you can write it in a way that relates it to your specific generative model for your data. That can tell you what features you ought to use if your model assumptions are correct. You can find the details in PRML section 4.2 (pdf). The essence is:

If your class labels are generated from a mixture of two classes with prior probabilities $\pi_1$ and $\pi_2$
If your covariates are generated from exponential-family distributions with densities proportional to $\exp(\beta_1^T\Phi(x))$ and $\exp(\beta_2^T\Phi(x))$

then the posterior probability that an observation $x$ is in class 1 is
$$\frac{\pi_1\exp(\beta_1^T\Phi(x)) }{\pi_1\exp(\beta_1^T\Phi(x)) + \pi_2\exp(\beta_2^T\Phi(x))$}$$. If you plug in known class labels and mangle this for a while with algebra, you end up with a logistic regression using features $\Phi(x)$. The intercept term of the logistic regression turns into a function of $\pi_1/\pi_2$.
If your data are Gaussian, then $\Phi(x)$ would generally include $x$ and $x^2$ and all the cross-terms if $x$ is multivariate. If there are special assumptions about equal covariance between classes, then the squared terms and cross terms may not be needed.


text:
Is there a way to tell which chromosome a gene is on, by looking at the "Chromosome/scaffold name"
I recently got a data set, from which I need to figure out which chromosome a gene is from, but the head of the data reads like:
Gene ID Description     Gene type       Gene End (bp)   Gene Start (bp) Strand  Associated Gene Name    Chromosome/scaffold name

while one entry of the data is like:
ENSG00000252303 RNA, U6 small nuclear 280, pseudogene [Source:HGNC Symbol;Acc:HGNC:47243]       snRNA   67546754        67546651        1       RNU6-280P       CHR_HG2128_PATCH

It looks like the last column is the information about which chromosome the gene comes from, but how to interpret the "CHR_HG2128_PATCH"?

What you are looking at is a patch, not an actual chromosome location. HG-2128 is an issue ID, something that the Genome Reference Consortium uses to track issue with references. Information regarding that specific patch is here. The issue is described as a "Possible misassembly or indel variation in GRCh38 within AL049860.8."
During the alignment step before you received your data, someone used a top level assembly version of the reference genome. This includes extra, accessioned scaffolds tacked on to the end of the file that contain the fix. These scaffolds represent what will be found in the next major release of the genome. The fixes aren't applied directly to the reference as soon as they are found due to the inevitable disruption of coordinates. So reads mapped to this FIX patch rather than the primary original scaffold in the primary assembly, on which the RNU6-280P gene resides.
The coordinates are slightly different. In the current GRCh38.p12
 assembly the RNU6-280P gene is located at chr6:67546651-67546754 (source). You can see that the coordinates are very close to those indicated in the table you provided (Gene Start, Gene End) but are indeed different.
Here is some more information on patches.


text:
Increase number of threads for GATK 4.0 HaplotypeCaller
I am using GATK version 4.0, I tried to use multiple threads for calling variants using HaplotypeCaller using following command
gatk --java-options -Xmx90G  -nt 28 HaplotypeCaller -I output.bam -R wheat_ref.fa -O final.vcf

and the error is
'-nt' is not a valid command

I have also used --num_threads but again the error was the same, as written here this parameter is available in previous versions of GATK but I don't know what's its status in current version of GATK i.e GATK 4.0. Whether it is deprecated or replaced by some other parameter?
Is there any way of calling multiple threads for variant calling with GATK 4.0?

I am assuming you are using the non-SPARK version of the method. To specify the number of threads you wish to use with HaplotypeCaller, include --native-pair-hmm-threads (documentation). This will only parallelize the pair hidden Markov models (pair HMM) process. This is an argument for the method HaplotypeCaller itself, so it should be located after the method is called as opposed to one of the --java-options.
I do recommend trying out the SPARK implementations of the GATK tools. The whole method will take advantage of the threads you specify. I've been running SPARK locally (i.e. not on Firecloud) on computational servers. I've been impressed with the enhanced performance.

I used this command now gatk HaplotypeCaller --native-pair-hmm-threads 28  -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf but it took same time and also I specified 28 threads and 28 was shown in requested threads but system monitor was showing only 3 ~ 4 threads are used as shown here https://imgur.com/a/JJJYmB7.


text:
Difference between samtools mark duplicates and samtools remove duplicates?
What is difference between samtools mark duplicates and remove duplicates ?  Is it necessary to mark duplicates before removing duplicates with samtools?

samtools rmdup and samtools markdup -r do the same thing. Without the -r flag samtools markdup only flags the duplicates.
You'll have to run samtools fixmate -m and sort the output to add ms and MC tags to prepare the file for markdup..
According to the current documentation, rmdup is obsolete. Please do not use rmdup.

TL;DR: just use markdup.
rmdup removes duplicates from BAM, while markdup, like Picard's MarkDuplicates, marks duplicates by default without hard removal – the latter is usually the desired behavior. In addition, markdup implements a better algorithm that takes care of more corner cases and gives more consistent results.
PS: to emphasize – I see no point to use samtools rmdup nowadays. It has been declared obsolete as @kohlkopf said.


text:
GRCm37-designed exon target enrichment, which reference to use?
I have exomes from 24 individual mice. The exomes are the product of Roche's SeqCap EZ HyperCap target enrichment kit. I see that the mouse exome design comes from mm9/NCBI37, the previous major reference genome released in July 2007.
Which reference should I use to align these reads? What artifacts may arise from mapping these reads designed from mm9 to the mm10 reference? I have already aligned the reads to mm10 and noticed an erroneous exon right in the middle of an intron, perhaps a target that was an exon in the mm9 release but revised out of the mm10 release. My goal is to investigate structural variation.

The only possible issue with using GRCm38 as the reference is that since GRCm37 was used in the bait design there are likely a few probes that will show artificial structural variations, due to errors in GRCm37. However, these will be present in all of your mice and if you use liftOver to look at the regions around the probes then such regions should become obvious. In fact, that'd be a convenient thing to do if you want to keep your (likely better) GRCm38 alignments, flank each probe site by a few hundred bases and then use liftOver to get the GRCm37 coordinates. If you go through the result you'll be able to find the false-positive structural variants due to the baits being designed with an older reference.
Alternatively, just use GRCm37 and hope that any issues with it don't result in false-negative SV calls.

I kept the GRCm38 aligned reads and used NCBI's remap tool on the target enrichment design alignment file from Roche. With the remapped target alignment file I was able to see that the exon in the middle of the GRCm38 enrichment was indeed a target! And good call on the false positive SVs being in all of the samples.

Can you clarify the SV issue? Regardless of how the probes are designed, wouldn't they only pull down the DNA fragments that physically exist?