text:
Strange per sequence GC content results
I would be grateful if someone could take a quick look at these FASTQC results. This is rna-seq paired-end data. From the FASTQC manual, an unusual distribution seems to be suggestive of contamination and a shift in the curve is suggestive of a systematic bias. GC content distribution both pre-alignment and post-alignment are strange.
Samples are paired end, strand specific and % of mapped reads is above 95% for all the samples. There is no adapter content also.
What could be the problem?

What are you trying to sequence? Something more specific than "poly-A amplified RNA" would be helpful to know if this is unexpected.

In the line of gringer. What specie are you sequencing ? Is there some background about GC distribution for this type of sequencing?  Did you check with other previous sequences to have such a problem? (It could be a termal microorganism and be totally normal or not...)

Its human species. except per sequence GC content and duplication levels of every sample other plots are totally fine. there is no adapter content also. samples are strand specific sequenced with poly-A selection. % of mapped reads is above 95% for all the samples. Paired end reads % average GC content is <50% or equal to 50%. % of duplicate reads is 65-70%

I see sequence bias in 5' and 3' end also in per base sequence plot. Is this GC content distribution due to bias selection?

What kit have you used to generate your libraries ?

Ribosomal rna depletion kit

Okay, so ribosomal depletion kits normally produce total RNA rather than polyA selected RNA. Normally you would do both ribo depletion and polyA selection. If you did ribo depletion rather than polyA selection, I would suggest you might be looking a contamination with ribosomal sequences. That would explain why you still see high mapping rates.

I'm not sure how normally people are combining ribo depletion and polyA, doing one obviates doing the other. However it's worth noting that ribo depletion is hit or miss regarding how well it actually works. I've seen cases with 1-2% rRNA afterward and cases with 20-30% rRNA, it depends a lot on how fresh the samples are.

@IanSudbery so is this common with ribo depletion kit? There is no adapter content also. Can I go further for with these samples?

Any kit can work inefficiently. You could check this by looking at how many reads map to ribosomal genes or ribosomal RNAs. Its not neccesasrily a problem, other than that you might end up with low counts in none ribosomal genes. I also be very careful when it came to normalisation in a DE experiment: large amounts of sequencing estate taken by a small number of genes can throw off normalisation for other genes.

Yes when I got the counts from the bam files after aligning to genome, I see most of the samples having approx. 25k-40k genes with 0 counts. Is this due to the reason what you said? How can I check reads mapping ribosomal genes or rRNAs?

@IanSudbery Hi, Could you please check my above comment. And could you please tell how to check reads mapping to ribosomal genes or rRNAs?

Sorry, our mapping pipeline does all this automatically and I had to dig out the code to remind myself how its done.
If you think you've got a very skewed counts distribution, i'd start by looking at the identities of the genes with the highest counts.
Alternatively to look directly at rRNA genes, you can download their annotations from UCSC or biomart with the correct biotypes. For the ribosomal proteins use go term GO:0003735. Then quantify with whatever you quantify with.
% reads mapping to protein coding genes is generally a pretty good QC metric.

Although generally I recommend people to be wary about the warnings that FastQC gives (it tends to be overly paranoid), the GC content graphs here do look odd. It's good that you've had a look at it with another program to confirm the observation. They should have a distribution that is fairly close to normal, a bit like this:

Without more information about your sample, I can't tell whether or not this is expected from your input. This looks like fairly substantial contamination, but could also be caused by biased selection during sample preparation (e.g. low quality sample, low yield, PCR amplification bias), or just looking at a multi-organism sample rather than a single-organism sample.
As a contamination check, I'd recommend running your RNASeq reads through a kmer-based classification tool (e.g. OneCodex, Centrifuge, BBSketch, or Sourmash). This should tell you the most likely place that contaminant sequences came from.


text:
Question about the dots on Quartile groups in boxplot
I have Microarray Normalized Expression data for a specific Gene. It looks like below in a dataframe B
SampleID    Gene    Type
Sample1     5.02    Tumor
Sample2     5.06    Tumor
Sample3     5.1     Tumor
Sample4     5.11    Tumor
Sample5    5.127    Normal
Sample6     5.12    Normal
Sample7    5.138    Normal
Sample8    5.149    Normal

I see that the minimum expression value in the table is 5.0 and maximum expression value is around 5.9.
I wanted to show the expression between two conditions with a boxplot and used following code.
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               add = "jitter", ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

This gave a plot like this
But when I remove the jitter from the code,
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

I see many black dots in the middle of the boxplot like this
May I know why I see those dots in the second boxplot after removing the jitter. Is there a way to avoid that?

I don't know from where does the ggboxplot function comes from or stat_compare_means? Also please, continue [edit]ing the question to clarify the question. It is still unclear to me if you want to remove some points, add them or something else

from the package ggpubr

I am confused by the wording in "May I know why I see those dots in the second boxplot after removing the jitter. Is there a way to avoid that?" Are you trying to ask  "Why did the points not disappear after I remove jitter?"

You add the points with geom_point(). Just remove it and you will get your "empty" boxplot.
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))

Unfortunately I couldn't use stat_compare_means(method = "t.test") and     stat_n_text(), but this two just add the labels.

yes, this is the mistake I'm doing. I removed geom_point(). It looks fine now. thanq

source: Stack Exchange

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

sort: votes


text:
*very* unbalanced group sizes for DE
I downloaded some publicly available RNA-seq data and want to compare those samples carrying a mutation (~4) against the rest (~800!). I ran both EdgeR and DESeq2, and the first results in an asymmetric volcano plot: skewed in one of the sides, meaning that it results in very large logFC that have large (NS) adjusted p-values. So next, in order to see if this asymmetry is persistant, I ran several random partitions of the data as groups (but keeping them highly unbalanced), obtaining the same behavior. Here an example of what I mean (I run it with 4, 20 and 100 to show how this effect decreases as I increase the small class size)

*I tried EdgeR-robust and I still get the same effect.
I do not see this effect when running DESeq2, only very few DE genes.  What method would you recommend using when comparing such different class sizes?

It'd be interesting to post this to the bioconductor site and see what Gordon Smyth says.

In general, you would expect roughly the same distribution on both sides of your volcano plot. The first two plots are concerning since there's so much on the right and not on the left. With such a small number of samples, it is hard to conclude if this is an artifact of not having enough samples, a technical artifact from the way the data was processed, or true biology. I remember reading that DESeq and other packages like a high number of samples for estimating dispersion or something, so this could also be contributing. I think you should find a better dataset (more samples) to prove or disprove your hypothesis.

I don't think this is a problem estimating dispersion: DEseq2 typically only requires a few replicates on either side (I've performed DE analysis on 3/3 with great results). Rather, the directionality could be due to the lack of power to detect significant differences in one direction vs the other.

@BenD. Wouldn't power be innately symmetric for cases like this? The power is finding a difference between two coefficient estimates (assumung a Wald test). Whether their difference is positive or negative would seem to be beside the point.

Random partitioning (as you have done) seems like a reasonable thing to do to work out what's going on. From the results you have suggested, it looks like DESeq2 is performing better than EdgeR, so use DESeq2. Can you show the volcano plots and the partition sizes?
Note: I'm a little bit biased in this, because DESeq2 is what I prefer to use (particularly where Salmon is not appropriate).

I think the same, I usually use DESeq2 too, but I was trying to repeat some analyses done by someone else with EdgeR. I edited my original post to include an image of the results.

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


text:
Script to run everything in a loop for extracting tar.gz files into fastq and to bam with alignment?
Original tar.gz files look like below:
UNCID_2668056.27d4c881-e906-437f-921d-467eb107798a.150108_UNC15-SN850_0402_AC5FTHACXX_2_GTTTCG.tar.gz
UNCID_2206882.58e94be5-2da2-4a96-944a-060f20216572.150108_UNC15-SN850_0402_AC5FTHACXX_3_CTTGTA.tar.gz
UNCID_2205121.2b3d596a-68ca-4c60-bb16-8c0cbd8fb190.140828_UNC11-SN627_0378_AC5FAWACXX_4_TAGCTT.tar.gz
UNCID_2209082.76a625e9-22b3-42b8-8813-27c133d0b248.110608_UNC11-SN627_0103_BD0CVVABXX_7.tar.gz

Before extracting the tar.gz files, I changed their names based on names.text file as below and using the command:
xargs -n2 mv -i <<<"$(sed '1d; s/$/.tar.gz/' names.txt)"

Tar.gz_filenames       Sample_Names
UNCID_2668056.27d4c881-e906-437f-921d-467eb107798a.150108_UNC15-SN850_0402_AC5FTHACXX_2_GTTTCG.tar.gz                Sample2
UNCID_2206882.58e94be5-2da2-4a96-944a-060f20216572.150108_UNC15-SN850_0402_AC5FTHACXX_3_CTTGTA.tar.gz                Sample1
UNCID_2205121.2b3d596a-68ca-4c60-bb16-8c0cbd8fb190.140828_UNC11-SN627_0378_AC5FAWACXX_4_TAGCTT.tar.gz                Sample3
UNCID_2209082.76a625e9-22b3-42b8-8813-27c133d0b248.110608_UNC11-SN627_0103_BD0CVVABXX_7.tar.gz                       Sample4

After changing the names, looks like below:
Sample1.tar.gz
Sample2.tar.gz
Sample3.tar.gz
Sample4.tar.gz

Using the tar.gz files after changing names I used the below code:
for i in /doc/*.tar.gz
do
    tar xvzf $i -C $TMPDIR
    for sample in $TMPDIR/*1.fastq
    do
        dir2="/destination/directory"
        base=$(basename $sample ".1.fastq")
        module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $TMPDIR/${base}.1.fastq -2 $TMPDIR/${base}.2.fastq | samtools view -Sb - > ${dir2}/${base}.bam
    done
done

When I run the above code, the output is like below:
150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_1.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_2.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_1.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_2.fastq
140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_2.fastq
140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_1.fastq
110608_UNC11-SN627_0103_BD0CVVABXX.7_1.fastq
110608_UNC11-SN627_0103_BD0CVVABXX.7_2.fastq
Warning: Could not open read file "/scratch/username/slurm-job.36320423/140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1
Warning: Could not open read file "/scratch/username/slurm-job.36320423/150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1
Warning: Could not open read file "/scratch/username/slurm-job.36320423/150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1

I dont mins whatever the name of fastq's in $TMPDIR but from $TMPDIR hisat2 should use those fastq's and give output bam files in /destination like below:
Sample1.bam
Sample2.bam
Sample3.bam
Sample4.bam

I would like to run a script doing tar.gz files -> fastq -> bam (everything in loop)

This should work for you:
for i in *.tar.gz
do
    dir2="/destination/directory"
    tar xvzf $i -C $TMPDIR
    files=($TMPDIR/*)
    module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done

(If this is the sort of thing you'll be doing often, it would be worth investing the time to learn a more robust approach such as Snakemake or Nextflow).
Edit:
If you want to run on all files at the same time, the easiest way is find / xargs
Change you script to this:
for i in $@
do
    dir2="/destination/directory"
    tar xvzf $i -C $TMPDIR
    files=($TMPDIR/*)
    module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done

Then run it like this:
find . -d 1 -name '*.tar.gz' | xargs -n 1 bash NAMEOFSCRIPT

This works for SGE. I'm not sure about the details of SLURM. You might have to figure out how to run each job in the background or use gnu parallel or something. As @Ian Sudbury says, pipeline software can take care of all these details for you.

Definately +1 for putting all this in some sort of pipeline. Many of these  will now handle automated distribution of many files to your SLURM cluster. Job tracking. Restarts after errors etc.

@heathobrien Hi, The above code uses one after another tar.gz file. It is taking huge time. How should I make the code work for all the tar.gz files at a time?

see updated answer

Even I haven't worked on SLURM before. Do you think if I give #SBATCH --array=1-6%3 will work? lets say I have 6 tar.gz files so 3 jobs at a time.

Your samples end with _1.fastq and _2.fastq, so you need to change this:
-1 $TMPDIR/${base}.1.fastq -2 $TMPDIR/${base}.2.fastq

to this:
-1 $TMPDIR/${base}_1.fastq -2 $TMPDIR/${base}_2.fastq

Thanq. ya this was the only error in my code.

Hi, The above code uses one after another tar.gz file. It is taking huge time. How should I make the code work for all the tar.gz files at a time?

Untar the files and use a workflow manager like snakemake. BTW, make sure you have sufficient resources to process many samples in parallel.

But as soon as Untar done for one file it should be the input for hisat2. Yes, I'm using 20 cpu cores, 1 week time and tmp 500G, memory usage 4G.

Then put the untar command into something like Snakemake too.

Never used a snakemake till now. May I know how I can give that?

Please google such things and post a new question if you have further questions.

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


text:
Detecting differentially expressed genes with foldchange >= 2 and FDR < 0.05
I'm using edgeR for differential analysis. Using glmTreat function I'm detecting differentially expressed genes between Tumor and Normal. I have set the arguments like below:
tr <- glmTreat(fit, contrast=contrast.matrix, lfc=2)
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

This gave me differentially expressed genes with logFC > 2 at FDR < 0.05 and genes with logFC < 2 at FDR < 0.05. [There are few genes with logFC 1.9, 1.89 also]
I have also tried with lfc=log2(2).
tr <- glmTreat(fit, contrast=contrast.matrix, lfc=log2(2))
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

To my surprise I see DE genes with logFC >=1, FDR < 0.05 and also logFC <=1 at FDR < 0.05.
Here is an example of DE genes I also found when I used lfc = log2(2)
GeneSymbols logFC   unshrunk.logFC  logCPM  PValue  FDR
DPYD-IT1    1.14493675  1.560389373 0.67468395  0.002972841 0.018704644
LINC00945   1.144288473 1.525191847 0.693845996 0.00651175  0.038064312
ITPKB-AS1   1.14046752  1.484991704 0.703887505 0.006620218 0.038585862
AC105206.1  1.122495266 1.4871905   0.699167758 0.006545299 0.038194633
LINC02066   1.120660287 1.558458881 0.66016102  0.004344743 0.026424298
AC090985.1  1.10235045  1.58968763  0.64343419  0.002745089 0.017426607
KCNJ6-AS1   1.100399688 1.717028637 0.613324799 0.002224039 0.01438451
LINC01494   1.090270138 1.490899364 0.676937984 0.006770608 0.039348027
AC099541.1  1.080957092 1.53817727  0.651356989 0.008449298 0.048356352
AL133396.2  1.012984207 1.512530231 0.63130643  0.008090259 0.046411928
AC011363.1  1.01207753  1.499708945 0.635740622 0.007226791 0.041797256

I'm very confused with these foldchange cutoffs.
Which one should I use to select differentially expressed genes based on fold change >=2 and FDR < 0.05?

Gordon Smyth has provided a comprehensive answer to your question below, but if any reader is interested in ranking genes by effect size check out TopConfects
(https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1674-7)

Fold-change >= 2 is the same as logFC (log2(fold-change)) >= 1, so your example is doing exactly what you want. logFC is generally easier to think in and work with than fold-change, since (aside from computational reasons) then increases and decreases in expression differ only in sign and are, thus, easier to compare in magnitude (e.g., it's easier to tell that logFCs of 1.5 and -1.5 are the same magnitude change than it is that 3 and 0.333 are the same magnitude change in the fold-change).
For prioritizing the results, typically you would then rank these filtered results by fold-change and relevance (either given background knowledge or for genes with known function).

Thanks a lot for the clear explanation. I guess in my example lfc = 2 is wrong analysis in this case if I need Fold-change >= 2.

The edgeR authors recommend that you use a relatively low logFC threshold for glmTreat such as lfc=log2(1.2). A lfc value as high as lfc=log2(2) is seldom required and only would only be appropriate for datasets with really large amounts of DE.
The first threshold you used of lfc=2 is equivalent to lfc=log2(4) and is too high for any dataset.
There's a few subtleties here:

First, there's the difference between logFC and FC as explained by Devon.

The lfc in glmTreat is not a fold-change cutoff! It is instead a testing threshold. All the significant genes will have actual logFCs well above the threshold.

The lfc threshold applies to the raw fold-changes (unshrunk.logFC in the output), not to the shrunk versions (logFC in the table).

The intent of glmTreat is to give a better gene ranking than would be provided by FDR or fold-change, so please don't re-sort the results by fold-change -- that would tend to defeat the purpose of the method.

I understand the documentation for glmTreat is a bit light, but you can get an idea of how to use it from our workflows.

https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

I understand the simple fold-change cutoff of FC > 2 with FDR < 0.05 is often used in the literature, but the edgeR authors do not recommend it and hence we don't provide an automatic way to get it from the edgeR pipeline.
For RNA-seq, fold-change cutoffs tend to prioritize genes with low expression and low counts.
The glmTreat approach is more subtle but we think gives better results, both in terms of biological relevance and statistical correctness.

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


text:
How to visualize genome track of gene in specific cell-lines?
I'm trying to make a plot showing genome tracks of specific genes in specific cell-lines of RNA-seq and Chip-seq data.
It should look something like this I have recently seen this encode, but in this I didn't find any MDAMB231 cell-lines.
May I know how to make such plots.
I got the Counts data from CCLE portal.
             BT-549 MDAMB231
RP11-554I8.2    0   246
LINC00707      50   107
   RBM17    17036   10819
  PFKFB3    18098   7668
RP11-393K12.2   0   0

Answer from @devon-ryan, converted from comments:
You'll need BAM or bigWig files for each of the tracks. If you can't find any you'll have to make them yourself from public fastq files. The fastq files can be obtained from ENA, SRA, or any of the other standard databases.

Like others have mentioned, these plots come from BAM, BED, and BigWig files.
Those files contain the data that you can give to tools like the UCSC Genome Browser, IGV,  HiGlass/Gosling, or other bioinformatics visualization software.
For the specific cell line you mention, you'll have to find the relevant RNA-seq and ChIP-seq data from some data repository, like CCLE, ENA, SRA, or GEO.
You may be to put in some work converting available files to one of these compatible formats before plotting them.

You can use screen.encodeproject.org
It collected a lot of data from multiple cell lines or tissues. And you can export tracks to UCSC browser or other browsers.


text:
Read counts from BAM file
I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:
samtools view -H IonXpress_input.bam

Output:
@HD VN:1.4  SO:coordinate   GO:none
@SQ SN:NR_039978    LN:984
@SQ SN:NM_130786    LN:1766
@SQ SN:NM_138932    LN:9293
@SQ SN:NM_033110    LN:2678
@SQ SN:NM_000014    LN:4678
@SQ SN:NM_144670    LN:5192
@SQ SN:NR_040112    LN:1201
@SQ SN:NM_017436    LN:2106
@SQ SN:NM_016161    LN:1771
@SQ SN:NM_015665    LN:1854

I guess these are from header. After these I see something like following:
@RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4

Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?

Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case).
For use in featureCounts my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).
If your data is properly aligned to hg19 your output should be something like this:
samtools view -H aligned2hg19.bam

@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10    LN:135534747
@SQ SN:chr11    LN:135006516
@SQ SN:chr12    LN:133851895
@SQ SN:chr13    LN:115169878
@SQ SN:chr14    LN:107349540
@SQ SN:chr15    LN:102531392
@SQ SN:chr16    LN:90354753
@SQ SN:chr17    LN:81195210
@SQ SN:chr18    LN:78077248
@SQ SN:chr19    LN:59128983
@SQ SN:chr20    LN:63025520
@SQ SN:chr21    LN:48129895
@SQ SN:chr22    LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566

So, I confirmed with the people who gave me these bam files. These are aligned to hg19 genome only.

Please show your output to these people.

But you mean the bam file I have doesn't have any chr, read count and other details? I was thinking the above posted one are from header and there is another command to get the readcounts and other details.

It probably has read counts, but not per chromosome but per @SQ sequence (which are RNA sequences). With samtools idxstats you'll see the read counts per @SQ.

Ok. I want to use featurecounts in R to get the read counts. I have the bam files and hg19.fasta and hg19.fasta.fai files. May I know how I should give the command to get the counts for all the samples. Thanq

Maybe you misunderstood, but you need to correctly align to hg19 first in order to do this in featureCounts. From your output it is clear that it is not aligned to hg19, but to a bunch of RNA sequences. Please discuss with your people that your data is not properly aligned to hg19 (show your output, this is proof), or try to align it yourself from fastq. When you have done that, you'll need a proper gtf annotation file from UCSC or ENSEMBL which is necessary for featureCounts. Good luck!

Your options are to either align to the genome (as @b.nota said) or filter using samtools however you like (e.g., for MAPQ) and then get counts with samtools idxstats.

When you are told they are aligned to hg19, I guess they mean tat are aligned to the hg19 version of the transcriptome, not hg19 itself. In this case you done need featureCounts to get read counts (and indeed, it would be hard/impossible to do so). samtools idxstats should give you want you are looking for.

Yes, thanks a lot for all. And yes it is aligned to hg19 transcriptome. And will give a try with idxstats.

One question. Do you I need to sort the bam files? samtools idxstats what are the options I should give in the command to get read counts?

Yes you'll need to sort bam files first. samtools idxstats aln.sorted.bam will give you the results, see here.

Yes, I saw the manual and tried like this and endup with an error.
samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam. This gave me [W::bam_hdr_read] EOF marker is absent. The input is probably truncated; [E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes;
[bam_sort_core] truncated file. Aborting.

The tmp dir has 6 input.sorted.nnnn.bam files.

Please start a new question for that error.


text:
sorting BAM file error using samtools
I have few bam files and would like to get read counts using
samtools idxstats

[Data is aligned to hg19 transcriptome]. To use that command I need a sorted bam file. So to sort them I gave the following command.
samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam

This ended up showing:
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
[E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes
[bam_sort_core] truncated file. Aborting.

And the tmp dir has 6 input.sorted.nnnn.bam files.

The cause of the error is that your input file is truncated or /tmp is running out of space. If you can do samtools view -H input.bam without error (reading the header also checks for the magic number at the end of the file) then it's like the that /tmp is filling up.

I did "samtools view -H input.bam" there is no error. It ended up showing all the header things which I asked here [https://bioinformatics.stackexchange.com/questions/3080/read-counts-from-bam-file/3081?noredirect=1#comment4889_3081]

Then either /tmp is full or otherwise one of the files there is somehow corrupt (you can test them with samtools view -H too).

Ok. looks like that file is corrupt. And I tried with other bam file and showed - [bam_sort_core] merging from 24 files... and I see a sorted.bam file. I used "samtools idxstats sorted.bam" - This is what I see - [bam_idxstats] fail to load the index.

Don't forget to samtools index sorted.bam

Ok. So, after the indexing I used idxstats and it showed 4 columns without any headers.                                                                                NR_039978 984 38 0
NM_130786 1766 1 0
NM_138932 9293 37 0
NM_033110 2678 2 0
NM_000014 4678 4 0
NM_144670 5192 18 0
NR_040112 1201 2 0
NM_017436 2106 76 0
NM_016161 1771 0 0

They are (in order): contig name, contig length, number of mapped reads, number of unmapped reads (this will typically be 0).

Ok. Thankyou. So, for DE analysis what I need is mapped reads. how do I get the gene names? Can I use edgeR for DE analysis using the mapped reads?

You mapped to the transcriptome, so they're the first column.

Yes, I see that. Is there any way to get the gene id/ transcript id from this results. And may I know the answer for DE analysis edgeR?

If you want to get annotation (like gene symbols) coupled to your accession numbers, please read the biomaRt user guide. If you have questions after reading it, please start a new question in stead of discussing it in the comments.


text:
Which module to select for Pathway analysis based on Module trait correlation and pvalue?
I have a total of 35 tumor samples classified into 4 subtypes. Subtype A, B, C, and D. I have RNAseq data. I'm interested in identifying modules related to each subtype with co-expression network analysis. I have read the WGCNA tutorial and create a trait.csv file like below:
Showing here some rows of the trait.csv file.

SamInfo <- read.csv("trait.csv")

Samples = rownames(datExpr);
traitRows = match(Samples, SamInfo$Samples);
datTraits = SamInfo[traitRows, -1];
rownames(datTraits) = SamInfo[traitRows, 1];

nGenes = ncol(datExpr);
#15648L

nSamples = nrow(datExpr);
#35L

#calculate eigengenes (1st principal component) of modules
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

#Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")",
                   sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));

#Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits),
               yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE,
               colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE,
               cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships"))

The Correlation between Subtypes and Modules Heatmap looks like below:

For functions identification (Pathway Analysis), from the above output for each subtype, should I select only the modules that were significantly positively correlated with subtypes, or Can I also use the significantly negatively correlated modules?
Is it good to use the modules that were not significant for pathway analysis?
In case if I can use positive and negative correlated modules, Can I combine those modules and use those genes for Pathway analysis?

It's a great question. Please consider the answer with a view to upvoting - or even accepting - the answer. It's good all round.

If you can run pathway analysis programmatically on all modules, yes, do it on all modules. Otherwise, select the modules with strongest associations with each subtype, both positive and negative. I would keep the modules separate because they may have very different biological functions.

Thanks for the answer. One small question. Lets say from the above image, for Subtype A I have selected modules salmon, grey60, tan (negative), and midnightblue (positive) for pathway analysis. I see that salmon module is also significantly positive with Subtype C. So, then the pathways I get from salmon module are related to which subtype?

Well, I don't see the figure anymore so I can only go by your description, but if salmon module is negatively correlated with A, it means (most of) the genes in the module are downregulated in subtype A and upregulated in (positively associated with) subtype C.

thanq. sorry for the heatmap with correlation values modules and trait information. Have some doubt whether here MEs0 = moduleEigengenes(datExpr, ?)$eigengenes I have to use dynamicColors which is before merging or moduleColors which is after merging. In the above image I used dynamicColors. 2) there are few modules like turquoise and blue with more than 8000 and 5000 genes. The number is huge. How do I reduce the number?

You can use either merged or un-merged colors; I would suggest using the merged ones.

Some modules are often very large which just reflects a factor (driving expression or co-expression) that strongly affects a lot of genes. This could be biological or technical, and could be important to you or a nuisance confounder. You could try to identify the factor e.g. from sample clustering (could be a batch effect) or by relating the MEs of largest modules to known sample characteristics (age, sex, race etc). If you identify a match, adjust the data for it, or simply adjust the data for the ME and run WGCNA again.

thanq very much. one more question. To visualize all the modules as network plot....Is it only cytoscape an option? Is there any other software or code to create network plot....


text:
Reads mapped to exonic, intronic and intergenic regions
After the alignment step I checked the rnaseq metrics of all the samples. Among 40 samples three samples show high percentage of reads mapped to intronic regions. What could be the reason?
Samples   Exonic            Intronic    Intergenic
Sample1 545479 (12.8%)  2512309 (58.8%) 1217201 (28.5%)
Sample2 372836 (8.3%)   2556934 (56.8%) 1573032 (34%)
Sample3 529618 (7.3%)   3934615 (54.5%) 2750720 (38.1%)

I also see that 35-40k genes among 58k genes having "zero" read counts. Is this due to contamination? What could be the reason? What does reads mapping to exonic, intronic, intergenic tell?
Update: I used hisat2 for alignment. I used human genome "grch38_snp_tran" from hisat2 website. Libraries were generated using ribosome depletion kit.

How did the overall alignment rate of those samples compare to the others that had more reasonable intronic alignment rates?

I see that all the 37 samples have overall alignment rate between 90-96% and for those 3 samples it is 69%, 79%, 70% overall alignment rate.

Sounds like you have degraded samples

U mean those 3 samples are degraded? What I should do now? Should I remove them for further analysis?

It’s a likely cause, but regardless you should exclude them.

Sure. the high percentage of mapping to intronic region is due to degradation?

Could you add the quality metrics of these files and other of the same run to the question?

There is the potential that these samples had low quality RNA or genomic DNA contamination (depending on the sample type and RNA isolation method)

What software did you get this output from?

As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.

I would expect at least 30% of reads from a total-cell, ribo-depleted RNA-seq to be exonic. Less suggests something when wrong.
As well as degradation, another explanation would be contamination with genomic DNA.
RNA obtained from the nuclear or cytoplasmic fraction of the cell might have a exonic content lower than 30%, while poly-A selected RNA would be expected to have higher, more like 66%.


text:
Mark BWA-SW split alignments in output for long reads
Is there a way to remove split alignments of single-end long reads from BWA-SW output? BWA-MEM has an option to include or flag split alignments in the result, but it seems BWA-SW method include them in the result unless the split occur in paired-end reads.

From BWA manual, it is advised to use BWA-MEM in your longs read case (unless you have a good reason not to):

BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate.

I guess that BWA-SW returns several lines per transcript (based on local alignments). In that case, the issue has already been reported. The user asks:

How can I differentiate between a read that has multiple hits in the reference and a read that has a split alignment?

The author replies:

In principle, it is hard to tell. Bwa-mem makes an arbitrary choice.
  Mappings flagged with 0x100 are multiple hits. Mappings flagged with
  0x800 are from a chimeric alignment.

So you could actually filter by flag. But I am afraid the output could include something you do not necessarily want to (since there is no specific flag for the split reads).
Maybe a better solution: you could switch to some mapper that will report the proper split in the hit of the read map, such as minimap2 from the same author. It is suitable for long reads and comes with a splice -x feature.
Then, the laziest solution I see here would be to remove the split alignments exploiting the sam format and in particular the CIGAR string. Meaning, to remove the sam fields containing a N in their CIGAR (aka, skipped region from the reference). Something like this would work:
cat alignment.sam | awk '$6 !~ /N/ {print $0}'

Hope this can be somewhat helpful.


text:
Meaning of BWA-MEM MAPQ scores
Does anyone know what the MAPQ values produced by BWA-MEM mean?
I'm looking for something similar to what Keith Bradnam discovered for Tophat v 1.4.1, where he realized that:

0 = maps to 5 or more locations
1 = maps to 3-4 locations
3 = maps to
  2 locations
255 = unique mapping

I'm familiar with the notion that MAPQ should be theoretically be related to the probability of an "incorrect" alignment (10^(-MAPQ/10)), although this is vaguely-specified enough that aligners in practice tend to actually just use something like the above.
Note that this is different than the BWA MAPQ scoring interpretation, because BWA-MEM gives MAPQ scores in the range $[0,60]$, rather than $[0,37]$ as has been established for BWA.

First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway.
Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the phred scaled probability of the alignment/base being wrong. It often amuses me that we question mapQ but take baseQ for granted. BaseQ is also scaled and discretized differently; even fewer people know how Illumina/pacbio/nanopore/historical sequencers estimates baseQ.
On the second question, Section 2 of the MAQ supplementary explains the theoretical aspects of mapQ, which is still correct today. Briefly, mapping quality consists of three components: 1) the probability of contamination; 2) the effect of mapping heuristics and 3) the error due to the repetitiveness of the reference. Only 3) can be modeled theoretically.
In case of bwa-mem, if we assume the matching score is 1, type-3 error is estimated with:
$$10/\log10\cdot[\log4\cdot(S_1-S_2)-\log n_{\rm sub}]$$
where $S_1$ is the best alignment score, $S_2$ is the second best and $n_{\rm sub}$ is the number of suboptimal alignments. Factor $\log 4$ comes from the scoring matrix. Factor $10/\log10$ is the Phred scale. This equation assumes gap-free alignment and is very close to Section 2.5.2 in the MAQ supplementary. It is ok-ish for short reads, but often overestimates for long reads. I am not aware of a practical approach in general cases. In addition to this method, you can estimate mapQ by read simulation: just try to find a function such that it fits empirical mapQ. Some have tried machine learning, too.

This has an ambiguous definition, at least because it doesn't define what happens when there are two equally similar and optimal matches in different parts of the reference index. If "second best" is allowed to have the same alignment score as "best", then this score reduces to $$10/\log10\cdot[-\log n_{\rm sub}]$$ in such a case.

Beautiful. Do you have a reference for the type-3 error estimation expression for bwa-mem?

@gringer Read the MAQ supplementary. It is an approximate on the condition that $S_1\gg S_2$. Furthermore, you do need to consider type-2 errors and pairing. The real estimate is full of ad hoc. ijoseph, that would the MAQ paper.

Those numbers are not arbitrarily picked (well... maybe 255/60/40 is arbitrarily picked).
To convert from log10 Q values like these (also used for error rates in FASTQ files) to probabilities, divide the number by 10, negate it, then raise 10 to the power of the result.
Another way of looking at it is to consider the decade to mean the number of 9s in the [mapping] accuracy, e.g.:
Q10: 90%
Q20: 99%
Q30: 99.9%
Q40: 99.99%
With the formula adjusted to make all the numbers in between a smooth transition.
With this transformation, a MAPQ score of 3 corresponds to a probability of an incorrect mapping of almost exactly 50%, which is as confident as you can be about a specific mapping if you know that it maps equally well to two locations:
> 10^-0.3
[1] 0.5011872

A MAPQ score of 2 is an incorrect mapping probability of 63% (i.e. close to 3 different locations), and a MAPQ of 1 is an incorrect mapping probability of 79% (close to 5 different locations). These fit with the numbers that you have described.
However, as you point out, BWA doesn't do this; it seems to ignore the number of mappable locations when considering the MAPQ score. I explore the reason why this is a problem in my question about multi-mapped reads.
I notice that there's a bit of an attempt to answer this in a similar previous question:

$C * (s_1 - s_2) / s_1,$
where $s_1$ and $s_2$ denoted the alignment scores of two best alignments and C was some constant.

See the answer by user172818 for more insight into the details of BWA-MEM MAPQ scores.

"Similar" question is generous, that's basically a duplicate of my question that I didn't do due diligence to find!

So it basically comes down to understanding this uncommented C code. Maybe I'll try contacting Heng Li directly. https://github.com/lh3/bwa/blob/master/bwamem.c#L952-L976


text:
Simulate fast5 MinION run
I'm planning on using Oxford Nanopore MinION sequencer to test a script that I'm working on. However, I would not like to run multiple MinION runs to validate the tool. I'd like to simulate various minION runs and get raw reads and all pore events in a fast5 file.  Do you know of any Oxford Nanopore (MinION) sequencing simulator that allows to get fast5 as output?
I looked into NanoSim but it seems like there is no option to get fast5 as an output.

Possible duplicate of Tools for simulating Oxford Nanopore reads

Hi Andresito. This question seems to be a duplicated of the linked question by gringer. However you could [edit] it to ask about the fast5 output if you would like to have an answer (and to keep open the question).

Using a simulator for example MinION output doesn't make sense when there are an abundance of freely-available FASTQ and FAST5 sequences for general purpose use (and especially for the purposes of testing out algorithms). Simulation of sequences can only ever express the components of DNA that are explicitly modelled, which is never going to fully represent reality.
As long as you're doing your own private research, you can use the FAST5 sequences from Clive Brown's own genome:
https://github.com/nanoporetech/ONT-HG1
Bear in mind the license attached to that data:
https://github.com/nanoporetech/ONT-HG1/blob/master/LICENSE.md

License to use for research use only (i.e. non-commercial) governed by the law of England and Wales.
You agree to this License by downloading the data on this site whether in whole or in part.
You cannot re-distribute modified copies of these data, or derivative works, without my express written permission. Where such permission is granted, you must cite the original work and its Copyrights and you must include this license with any work that includes any of my genome.

For a more open license, you can use the nanopore-wgs reads here (derived from human cell line NA12878):
https://github.com/nanopore-wgs-consortium/NA12878

We encourage the reuse of this data in your own analysis and publications which is released under the Creative Commons CC-BY license. Therefore we would be grateful if you would cite the reference below if you do.

If you don't want human data, there's mouse parasite data here:
http://www.ebi.ac.uk/ena/data/view/PRJEB20824
And E. coli data here:
https://www.ebi.ac.uk/ena/data/view/PRJEB22772
... and if you just want 8h of sensor noise as a FAST5 file, I've got that here in one of my DropBox folders (as a ~70MB file):
https://www.dropbox.com/s/20j1anlz0motrs4/MINICOL233_20170118_FN__MN15625_sequencing_throughput_ONLL3127_84599_ch320_read10302_strand.fast5?dl=0
Reminder: if you do ever present the results of analyses based on a data source, it's a good idea to acknowledge the source (even if not expected). If nothing else, it will let the creator know that someone cares about their data, and perhaps encourage them to release more in the future.

Thank you very much @gringer. I will revise the provided sources and acknowledge accordingly. You are great!

I wrote a tool that will rerun an experiment. So all you would need to do is just download an ONT dataset from ENA or something and you can just rerun it as many times as you want.
https://github.com/mbhall88/taeper


text:
Can we use non base-called fast5 files in poretools?
I run the MinION MinKNOW without the live base-calling option. We know there is Metrichor and Albacore to perform base-calling after this process. However, I have not done any base-calling yet.
My question is: Is it possible to use directly the fast5 files with poretools to extract a fastq file without doing any previous base-calling process? I tried it and I get an empty fastq file.
My reads directory only contains the fast5 folder and I run:
poretools fastq fast5/ > output.fastq
Any ideas why do I get an empty file? What is the difference between base-called fast5 and non-base-called fast5 files?
Now I am trying to do base-calling with Albacore to see if I get a fastq file.

You’ve edited the title - I should add that poretools can work with non base-called FAST5 files with some of its options, but not the fastq option. Some of the metadata options will function.

No, poretools does not do basecalling. The poretools fastq command can be used to extract the FASTQ information from the basecalled FAST5 file (via MinKNOW live-basecalling or albacore). Alternatively, both of these basecallers can export a FASTQ file directly if desired.
The difference between a basecalled FAST5 and a nonbasecalled FAST5 is the presence of the FASTQ within the file. This can be directly viewed using an HDF5 viewer, like h5dump from the commandline, or HDFView from the GUI.

Albacore can be used to create a fastq, directly from your fast5 files. As such there is no more reason to use poretools, expect perhaps for getting some metadata. With regards to metadata I would like to add that the sequencing_summary.txt generated by albacore also contains a lot of metadata.


text:
cellranger-arc throws an error when demultiplexing
I have a multiome dataset from a 10x run I'm trying to demultiplex using cellranger-arc mkfastq.
The demultiplexing uses the cellranger-arc mkfastq command with a samplesheet in either of two forms, a simple one using a csv file, or a more extensive one called IEM (Illumina Experiment Manager).
Both formats throws the same error, so I don't think it is the file itself.
The command I'm using to start the demultiplexing is
cellranger-arc mkfastq --id=P786_GEX --run=/fs/pool/pool-bcfngs/Raw_data/P786_240424_NB500982_0590_AHW55LAFX5/ --samplesheet=/fs/pool/pool-cox-projects-bioinformatics/MPIB_NGS/P786/P786-gex-samplesheet-1.0.0.csv --use-bases-mask=Y28n*,I10,Y10n*,Y82 --output-dir /fs/pool/pool-bcfngs/Raw_data/P786_GEX/

When using the simple csv file the only difference is one has to use the parameter --csv instead of --samplesheet.
The samplesheet itself is a comma-seaprated value file and looks like that:
Lane,Sample,Index
*,RNAseq_A,SI-TT-B1
*,RNAseq_R,SI-TT-C1

The * means that the sample can be found in all four lanes. This was run on a NextSeq500 machine, which can't separate the four lanes.
The error I get is this is as such:
$ cellranger-arc mkfastq --id=P786_GEX --run=/fs/pool/pool-bcfngs/Raw_data/P786_240424_NB500982_0590_AHW55LAFX5/ --samplesheet=/fs/pool/pool-cox-projects-bioinformatics/MPIB_NGS/P786/P786-gex-samplesheet-1.0.0.csv --use-bases-mask=Y28n*,I10,Y10n*,Y82 --output-dir /fs/pool/pool-bcfngs/Raw_data/P786_GEX/
cellranger-arc mkfastq (cellranger-arc-2.0.2)
Copyright (c) 2021 10x Genomics, Inc.  All rights reserved.
-------------------------------------------------------------------------------

Auto-detecting ATAC/GEX: ATAC mode - I2 is 16 bases
Martian Runtime - v4.0.5
Serving UI at http://hpcl9301:42121?auth=cw0ukK_kaNpv_kf2sTe_pJRF0OXf3B7gaVweU23Pjh0

Running preflight checks (please wait)...
Checking run folder...
Checking RunInfo.xml...
Checking system environment...
Emitting run information...
Checking read specification...
Checking samplesheet specs...
2024-04-29 13:27:53 [runtime] (ready)           ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.PREPARE_SAMPLESHEET
2024-04-29 13:27:53 [runtime] (run:local)       ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.PREPARE_SAMPLESHEET.fork0.chnk0.main
2024-04-29 13:27:54 [runtime] (chunks_complete) ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.PREPARE_SAMPLESHEET
2024-04-29 13:27:54 [runtime] (ready)           ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET
2024-04-29 13:27:54 [runtime] (run:local)       ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET.fork0.split
2024-04-29 13:27:54 [runtime] (split_complete)  ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET
2024-04-29 13:27:54 [runtime] (run:local)       ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET.fork0.chnk0.main
2024-04-29 13:32:16 [runtime] (chunks_complete) ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET
2024-04-29 13:32:16 [runtime] (run:local)       ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET.fork0.join
2024-04-29 13:32:16 [runtime] (join_complete)   ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.BCL2FASTQ_WITH_SAMPLESHEET
2024-04-29 13:32:16 [runtime] (ready)           ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.MERGE_FASTQS_BY_LANE_SAMPLE
2024-04-29 13:32:16 [runtime] (run:local)       ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.MERGE_FASTQS_BY_LANE_SAMPLE.fork0.split
2024-04-29 13:32:17 [runtime] (failed)          ID.P786_GEX.MAKE_FASTQS_CS.MAKE_FASTQS.MERGE_FASTQS_BY_LANE_SAMPLE

[error] Pipestance failed. Error log at:
P786_GEX/MAKE_FASTQS_CS/MAKE_FASTQS/MERGE_FASTQS_BY_LANE_SAMPLE/fork0/split-u7b0e2f8540/_errors

Log message:
Traceback (most recent call last):
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/external/martian/adapters/python/martian_shell.py", line 659, in _main
    stage.main()
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/external/martian/adapters/python/martian_shell.py", line 610, in main
    lambda: self._record_result(lambda: self._module.split(args))
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/external/martian/adapters/python/martian_shell.py", line 589, in _run
    cmd()
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/external/martian/adapters/python/martian_shell.py", line 610, in <lambda>
    lambda: self._record_result(lambda: self._module.split(args))
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/external/martian/adapters/python/martian_shell.py", line 593, in _record_result
    self._result = cmd()
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/external/martian/adapters/python/martian_shell.py", line 610, in <lambda>
    lambda: self._record_result(lambda: self._module.split(args))
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/mro/tenkit/stages/make_fastqs/merge_fastqs_by_lane_sample/__init__.py", line 44, in split
    iem_df = tk_sheet.file_get_iem_data_frame(args.samplesheet_path)
  File "/fs/gpfs41/lv07/fileset07/home/b_cox/yeroslaviz/cellranger/cellranger-arc-2.0.2/lib/python/tenkit/samplesheet.py", line 351, in file_get_iem_data_frame
    raise ValueError("Invalid IEM samplesheet format: %s" % path)
ValueError: Invalid IEM samplesheet format: /fs/pool/pool-bcfngs/Raw_data/P786_GEX/MAKE_FASTQS_CS/MAKE_FASTQS/PREPARE_SAMPLESHEET/fork0/chnk0-u7b0e2f8439/files/samplesheet.csv

I don't know wehre to look for the error. I would appreciate any help in this regards.
EDIT
After swbarnes2 comment I added this here:
so, yes the samplesheet is correct. There are two different kinds of samplesheets and Illumina and 10x both recommend using the simple format of a csv file. But I did try though both of them.
I also tried this one here:
[Data],,,,,,,,,,
Lane,Sample_ID,Sample_Name,Sample_Plate,Sample_Well,I7_Index_ID,index,I5_Index_ID,index2,Sample_Project,Description
1,RNAseq_A,RNAseq_A,,,SI-TT-B1,SI-TT-B1,SI-TT-B1,SI-TT-B1,P786,
2,RNAseq_A,RNAseq_A,,,SI-TT-B1,SI-TT-B1,SI-TT-B1,SI-TT-B1,P786,
3,RNAseq_A,RNAseq_A,,,SI-TT-B1,SI-TT-B1,SI-TT-B1,SI-TT-B1,P786,
4,RNAseq_A,RNAseq_A,,,SI-TT-B1,SI-TT-B1,SI-TT-B1,SI-TT-B1,P786,
1,RNAseq_R,RNAseq_R,,,SI-TT-C1,SI-TT-C1,SI-TT-C1,SI-TT-C1,P786,
2,RNAseq_R,RNAseq_R,,,SI-TT-C1,SI-TT-C1,SI-TT-C1,SI-TT-C1,P786,
3,RNAseq_R,RNAseq_R,,,SI-TT-C1,SI-TT-C1,SI-TT-C1,SI-TT-C1,P786,
4,RNAseq_R,RNAseq_R,,,SI-TT-C1,SI-TT-C1,SI-TT-C1,SI-TT-C1,P786,

It gave me the same error
more can be found here or here:

I have found the solution to my problem. The only thing missing from the command I used was the instruction to regard the run as a ATAC/GEX step. This is done by adding the parameter --filter-dual-index for GEX to the command.
The complete, working commands for both parts of the same data set are therefore these:
cellranger-arc mkfastq --id=ATAC \
                       --run=P786_AHW55LAFX5/ \
                       --csv=ATAC.csv \
                       --filter-single-index \
                       --use-bases-mask=Y50,I8n*,Y16,Y82 \
                       --output-dir ATAC/

cellranger-arc mkfastq --id=GEX \
                       --run=P786_AHW55LAFX5/ \
                       --csv=GEX.csv \
                       --filter-dual-index \
                       --use-bases-mask=Y28n*,I10,Y10n*,Y82 \
                       --output-dir GEX/

the number in the --use-bases-mask depends on the run parameters (see RunInfo.xml file) and must be adjusted accordingly.

I don't think that samplesheet is right.  Just make 4 entries per sample, one per lane.

I disagree. At least I don't see anything wrong with the format. Pls see my edit to my question.


text:
How to determine the primary Uniprot accession number from a set of accession numbers?
Given a list of Uniprot IDs that are linked to an Ensembl gene ID, is there a way to systematically determine the UniProt accession corresponding to the protein product of the gene's primary transcript with no other information?
According to ExPasy

Researchers who wish to cite entries in their publications should
  always cite the first accession number. This is commonly referred to
  as the 'primary accession number'. 'Secondary accession numbers' are
  sorted alphanumerically.

But what if the order has been jumbled or compiled from a different source and resorted?
For example:
Ensembl:
ENSMUSG00000035642
Uniprot:
Q8R0P4, Q8CF11, D6RJK8, D6RJJ4, D3Z442, D3Z1Q3, D3YZD8, D3YY39, D3YX09, D3YWY5
This question is cross listed on the Biology stack-exchange.

If I understand this correctly, each Uniprot entry has the same primary accession number (which you can find in the "Entry Information" tab of the Uniprot entry; ctrl+f "primary"). Some of the IDs that you give here describe different sequences (which thus have different Uniprot entries and different primary accession numbers). I have checked the first 5 IDs you show here. They all map to different sequences, and should have different accession numbers. They are isoforms, mutated or truncated sequences and God knows what else, but they are definitely not the same protein sequence.

I never claimed these are representative of the same protein squence. The question above states: Given a list of Uniprot IDs that are linked to an Ensembl ID

No, you didn't :) What I am trying to say is that the way I understand it there shouldn't be an 1:1 mapping between Ensembl and UniProt in the first place. They describe different things. Am I missing something? It seems to me like you are trying to pick one protein ID to connect to one gene ID.

@galicae didn't do what? I'm trying to determine if there is a simple procedure to identify the primary Uniprot accession number from a list of associated accession numbers with no other outside information. I don't think I can state it more clearly.

@JamesDraper do you mean the primary UniProt accession of the protein product of the primary transcript associated with that gene? Or do you mean the primary UniProt accession of each of the protein isoforms encoded by the gene? Please [edit] your question and clarify what you need since there is often (usually) a one-to-many relationship between genes and proteins.

@JamesDraper you said you "never claimed these are representative". I agreed :)

@terdon the first option you mentioned. Feel free to edit the post as you see fit.

@galicae I'm glad that we are on the same page.

@JamesDraper OK, and how do you need to do this? Manually? Does it need to be scalable to hundreds of accessions? You mention "no other information", does that mean that going from the Ensembl Gene ID, to the Ensembl Transcript ID for the primary transcript and then getting the UniProt accession(s) associated with that is not an option?

@terdon honestly I was was hoping that this question had a simple obvious answer that I just didn't know like: sort in in descending order and take the first member of that list. Or if the something like the primary accession could be recognized by  some kind of special regex. However it has become abundantly clear that the answer to the question is no there is no way to determine the primary accession from the a list of associated accession numbers with no outside info. If anyone would like to restate that as answer I'll accept and upvote.

Reposting my answer from the related ticket in the Biology section:
I think there is an issue with the terminology. The "primary" accession number, is the first accession number in cases where an entry has more than one accession number, as described in http://www.uniprot.org/help/accession_numbers:

Entries can have more than one accession number. This can be due to
  two distinct mechanisms:
a) When two or more entries are merged, the accession numbers from all entries are kept. The first accession number is referred to as the

‘Primary (citable) accession number’, while the others are referred to
  as ‘Secondary accession numbers’. These are listed in alphanumerical
  order.
b) If an existing entry is split into two or more entries (‘demerged’), new ‘primary’ accession numbers are attributed to all

the split entries while all original accession numbers are retained as
  ‘secondary’ accession numbers.
Example: P29358 which has been ‘demerged’ into P68250 and P68251.

Both reviewed and unreviewed entries can have primary accession numbers.
What you probably mean, as previous posters understood, are accession numbers of reviewed entries as opposed to unreviewed ones.
In that case, you can indeed add "reviewed:yes" to your query, e.g. when you are using the UniProt ID mapping, http://www.uniprot.org/help/uploadlists

Consider the circumstances in which entries will have more than one accession (taken from the same page you linked in the original post).

Entries will have more than one accession number if they have been merged or split. For example, when two entries are merged into one, the accession numbers from both entries are stored in the AC line(s).
If an existing entry is split into two or more entries (a rare occurrence), the original accession numbers are retained in all the derived entries and a new primary accession number is added to all the entries.

As mentioned in the comments, this does not imply that there is or should be a single "primary" UniProt accession associated with (for example) an Ensembl gene ID. In some cases, a single protein record may have multiple accessions due to annotation changes, and the suggestion from your original post is to only report the primary accession. However, there may be multiple proteins associated with an Ensembl gene ID, and each protein will have its own accession (or set of accessions).
Choosing a single representative isoform for a gene is an entirely different question.

Not with no other information, no. There is no way of looking at a list of accessions and guessing which one is primary. Primary accessions are just accessions like any other; they don't have any specific characteristics.
This is further complicated by the fact that each gene can (and usually does, for higher eukaryotes) have multiple transcripts, each of which can produce a different protein product. So in order to go from gene identifier to a protein, you need to first find the relevant transcript.
So, as a first step, get the primary transcript IDs for each of your genes and then retrieve the protein accessions associated to them.


text:
QMean standalone evaluation
I want to use a script to evaluate a protein structure with QMEAN scoring.
I can't understand why the script suggested by SWISS-MODEL (access to QMEAN by json object) runs only if I put the commands on the shell. If I build a module, the result of the second request (request.get), return a QUEUEING status, not COMPLETE status.
Here is the code that I have been using:
import json
import requests

def qmean_page(structure):
""" return the Qmean4 score from SWISSMODEL QMEAN tool"""

    url = "https://swissmodel.expasy.org/qmean/submit/"

    # Request the resource
    response = requests.post(url=url, data={"email": "XXX@XXX.XX"} ,files={"structure": open(structure,'r'),},timeout=3,verify=True)

    json.dumps(response.json(), indent=4, sort_keys=True)
    download_json = str(response.json()['download_url'])

    return download_json

def qmean_score(download_url):
# Request the resource
    r = requests.get(download_url, verify=True)

    scores = json.dumps(r.json(), indent=4, sort_keys=True)
    print scores
    jdata = json.loads(scores)
    zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
    return zscore

s = 'P02489.B99990005.pdb'
url_d = qmean_page(s)
print url_d
score  = qmean_score(url_d)
print score

CASE A: If I comment the last 5 line, I can import qmean.py (the function file) in python terminal.
>>> import qmean
>>> s = 'P02489.B99990005.pdb'
>>> url_d = qmean.qmean_page(s)
>>> print url_d
https://swissmodel.expasy.org/qmean/project/TgEjha.json
>>> score  = qmean.qmean_score(url_d)
{
"input_data": {
    "info": {
        "model_001": {
            "chain_mapping": {
                " ": "A"
            }
        }
    },
    "meta": {
        "created": "2018-01-18T10:04:29.978",
        "email": "XXX@XXX.XX",
        "project_name": "Project",
        "seqres_uploaded": false
    },
    "method": "QMEANDisCo",
    "models": [
        {
            "error_status": null,
            "global_scores": {
                "acc_agreement": {
                    "norm": 0.5606936416184971,
                    "zscore": -1.0459143647477258
                },
                "all_atom": {
                    "norm": -0.012626634534455788,
                    "zscore": -2.1825427852787604
                },
                "cbeta": {
                    "norm": -0.004260555683274807,
                    "zscore": -1.9157144117276892
                },
                "qmean4": {
                    "norm": 0.6275190794418422,
                    "zscore": -3.5006904416158258
                },
                "qmean6": {
                    "norm": 0.6339071011931049,
                    "zscore": -3.0507104571131918
                },
                "solvation": {
                    "norm": -0.609099435253941,
                    "zscore": -1.2499924754236265
                },
                "ss_agreement": {
                    "norm": 0.6706229586994028,
                    "zscore": 0.05371252044245716
                },
                "torsion": {
                    "norm": 0.0033300496230023145,
                    "zscore": -2.7188811736675764
                }
            },
            "model_pdb": "https://swissmodel.expasy.org/qmean/project/TgEjha/model_001.pdb",
            "modelid": "model_001",
            "name": "P02489.B99990005.pdb",
            "seqres": [
                {
                    "atomseq": "MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSGISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRLPSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS",
                    "chain_name": " ",
                    "name": "seq_chain_0",
                    "sequence": "MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSGISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRLPSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS"
                }
            ]
        }
    ],
    "options": {
        "qmeanbrane": false,
        "qmeandisco": true
    },
    "project_type": "default",
    "sequences": []
},
"status": "COMPLETED"
}
>>>

CASE B: If I uncomment the last 5 lines and run 'python qmean.py', it doesn't work.
  andrea$ python qmean.py
  https://swissmodel.expasy.org/qmean/project/zrcLCK.json
  {
  "status": "QUEUEING"
  }
  Traceback (most recent call last):
  File "qmean.py", line 31, in <module>
  score  = qmean_score(url_d)
  File "qmean.py", line 24, in qmean_score
  zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
  KeyError: 'input_data'

If I print 'r' variable of 'qmean_score' function, It returns:
  <Response [200]>
  {
  "status": "QUEUEING"
  }

You need to use if __main__: if you want to use that file as a module and in command line. But what is your question? You seem to have solved your problems.

No, it is not solved because I could have a result like case A also when I call my function as case B

I resolve the problem by using a while loop.

maybe you could post your answer and accept it, for people trying to figure out what was the issue, or delete your question if you think that it won't be useful for other people.

So this is connected to how QMEAN is implemented as a server. It also could be true to many bioinformatics resources so I think this is important.
Posting to QMEAN "submit" triggers the job that is processed on our HPC cluster. That takes time and in the meantime JSON response just tells you that the job is QUEUEING/RUNNING or when it is ready it tells you COMPLETED. The answer suggested by you is the correct assuming you wish to wait for the job completion (which might take quite a lot).
So just check in while loop when the status is COMPLETED and fetch the result eg.
import time
import json
import requests
def qmean_score(download_url):
    # Request the resource
    zscore = None
    while True:
        r = requests.get(download_url, verify=True)

        scores = json.dumps(r.json(), indent=4, sort_keys=True)
        jdata = json.loads(scores)
        if jdata["status"] == "COMPLETED":
            zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
            break
        else:
            time.sleep(1000)
    return zscore

But this doesn’t differ between interactive and batch usage. OP just got lucky when running the code from an interactive session.

It does in a sense that one can submit many targets, save them in eg. file and after a while, one can return to them and retrieve all the date with JSON. It is just a convenience API for people that would like to do it in the command line and not doing this by hand. It is also used in CASP and CAMEO.


text:
Primer-Blast Custom Database [error]
I'm using Primer-Blast to compare a PCR template against a custom database of reference sequences, using Primer-Blast's Custom database functionality.
However whenever I try to upload a FASTA file as the custom database and then run Primer-Blast, I get the following error:
Exception error: Fasta Reader: sequence id ends with 80 valid nucleotide characters. Was the sequence accidentally placed in the definition line?
No matter how I try to format my FASTA file, I get this error.
Does anyone know how to properly format a FASTA file to have it accepted as a custom database of reference sequences to compare against the PCR template? Or is this a bug in the software?
Here's an example FASTA file, to upload both for the template and the custom database.
>test1
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTGC

Could you provide a sample of your FASTA file, pls?

@mrhd just added

@mrhd It's almost like it's not recognizing the newline escape character at the end of the sequence ID.

That might actually be the cause of the error. Did you create the FASTA on Windows or Linux/Mac?

@mrhd The actual FASTAs I want to use were downloaded from GISAID. However I am on a Mac.

Ok so it looks like the workaround for this is to copy-paste the template into the text box, and then upload your custom database. If you upload both the template and the custom database, you will get the error.

I just copied your FASTA sequence into a file (on Mac OS) and uploaded it to Primer BLAST, the search did not give an error.
If the error is caused by the line endings, maybe you can fix your downloaded files with a tool like dos2unix.
Update:
On a whim, I tried adding a description to the FASTA header:
>test1 description
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTG

This time, there was no error.
Does that work for you as well?

You uploaded it both as the PCR template and as the custom database? Because it works fine when you just upload it as the template, and then set the database to something like Refseq representative genomes. It's only when you set the database to custom that it errors out.

Sorry, you're right. When uploading the sequence as database, I get the same error. -.-

The error does not appear when pasting the sequence into the text box under "Custom", though. Strange...

Ok, for some crazy reason, I'm able to get it to work when I copy-paste the template but upload the custom database. But it doesn't work when you upload both the template and the custom database, even though it does work if you upload the template but don't use the custom database option.

Very odd indeed..

I did it by copy and paste. Here's the top line of the output:
    Primer pair 1
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  TGTCGTTGACAGGACACGAG    Plus    20  148 167 59.97   55.00   5.00    2.00
Reverse primer  TTACCTTTCGGTCACACCCG    Minus   20  264 245 59.97   55.00   3.00    2.00
Product length  117
Primer pair 2
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  ATCTGTGTGGCTGTCACTCG    Plus    20  82  101 60.04   55.00   5.00    3.00
Reverse primer  ACTCGTGTCCTGTCAACGAC    Minus   20  168 149 59.97   55.00   4.00    3.00
Product length  87
Primer pair 3
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  TACTGTCGTTGACAGGACACG   Plus    21  145 165 60.00   52.38   7.00    3.00
Reverse primer  ACACCCGGACGAAACCTAGA    Minus   20  251 232 60.54   55.00   4.00    2.00
Product length  107
Primer pair 4
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  GCTTACGGTTTCGTCCGTGT    Plus    20  192 211 60.94   55.00   4.00    0.00
Reverse primer  GTTGAAACCAGGGACAAGGC    Minus   20  292 273 59.33   55.00   4.00    2.00
Product length  101

A command line solution would work, based on the bug reported by the OP. However there's nothing in NCBIs e-utilities and Biopython doesn't do Primer-Blast, but will do Blast. You could do it Perl via the LWP Agent, e.g. here.


text:
BLASTn query coverage discrepancy
Using the BLASTn webtool with,
Query Seq: CAATTGCGAAGCCACATTATT
Subject Seq: CP031332.1
BLASTn returns a single hit to CP031332.1 with query coverage at 100%.
However if one looks at the alignment, the query coverage is clearly 17/21 ≈ 81%. Could someone explain to me the reason for this discrepancy?
Edit,
Web version: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome
Output: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=RCK3F55Y114

To expand on the excellent answer by haci:
BLAST reports (a) the summary of alignments and (b) individual alignments, see below. Note that your link to the results expired (as NCBI BLAST results links tend to do). It is considered best practice to post the text version of the relevant portion of the results in the question itself.
The summary provides an overview of all the query sequence individual hits in the subject sequence. Because the individual hits (alignments, or ranges) together cover 100% of the query and each of these has the max identity = 100%, that is what it reports in the summary. It is not very intuitive for many users.
You should always look at the individual hits in detail. The above 100% coverage in the summary consists of ranges 1 and 14, shown below. The top hit (range 1) covers only 17 out of 21 nt with 100% identity, and this is the most relevant information from BLAST. The other hit (range 14) covers 10 out of 21 nt with 100% identity. This is the cause of the discrepancy.
These 2 alignments are not close together in the subject sequence. Depending on your goal, choose one or more of the individual alignments and consider alignment length and % identity of those hits individually, rather than summarized together as in the BLAST report summary. For most applications, the top hit (range 1) would be the most relevant alignment.
Program: BLASTN
Subject:Neisseria meningitidis strain M22814 chromosome, complete genome ID: CP031332.1(nucleic acid) Length: 2190201
Query #1: seq1 Query ID: lcl|Query_56011 Length: 21

Sequences producing significant alignments:
                                                                  Max    Total Query   E   Per.
Description                                                       Score  Score cover Value Ident  Accession
Neisseria meningitidis strain M22814 chromosome, complete genome  34.2   57729 100%  0.001 100.00 CP031332.1

Alignments:
>Neisseria meningitidis strain M22814 chromosome, complete genome
Sequence ID: CP031332.1 Length: 2190201
Range 1: 2074883 to 2074899

Score:34.2 bits(17), Expect:0.001,
Identities:17/17(100%),  Gaps:0/17(0%), Strand: Plus/Plus

Query  1        CAATTGCGAAGCCACAT  17
                |||||||||||||||||
Sbjct  2074883  CAATTGCGAAGCCACAT  2074899

...

Range 14: 354292 to 354301

Score:20.3 bits(10), Expect:15,
Identities:10/10(100%),  Gaps:0/10(0%), Strand: Plus/Minus

Query  12      CCACATTATT  21
               ||||||||||
Sbjct  354301  CCACATTATT  354292

Your blast result is not composed of a single alignment of all of the length of your query vs the subject but from a series of alignments: in each a subsequence of your query is used and these alignments are sorted by "E value", "score", "percent identity", ...
Now in your case, you are talking about the first hit, in which a subsequence of your query (positions 1-17) match perfectly to a portion of the subject. Query coverage is not calculated from this first hit though, it takes all of the alignments into account and again in your case all 21 bases were used in one or more alignments. I believe it would be fair to filter out all but the first hit given the E-values and re-calculate the query coverage as 17/21.


text:
blastn returning an inferior alignment
In my work I've found that in circumstances where there is a mismatch near the edges of the query sequence, blastn prefers to return a shorter contiguous alignment, rather than allowing for a mismatch near the 5' or 3' end of the query sequence.
For example, consider the subject sequence (bold indicates alignment with mismatch in brackets):
GACCAACCTGCTGGATCTGCTGATCGGCG[G]TTTGCGAAATA
and query sequence:
CTGGATCTGCTGATCGGCG[A]T
Blasting this will return the 19bp contiguous alignment upstream of the mismatch, rather than returning a 21bp alignment with a single mismatch at the G/A base at position 20 in the query.
I'm using the command line version of blastn, and I've tried messing with the penalties as well as setting -task blastn-short, but nothing seems to fix this issue.
I speculate that the various gap and mismatch penalties are working out such that a shorter contiguous alignment is preferred over a slightly longer alignment with a mismatch, but I haven't been able to figure out how to modify the penalties in such a way as to change this preference.

Nope, altering the weights is the solution as you said. So please do answer your own question —this a totally encouraged behaviour.

Here's my working solution,
So I found a choice of parameters that resolves this, -reward 3 -gapopen 5 -gapextend 5.


text:
Software recommendation: find DNA sequence distribution over entire transcript
I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms.
I have already found the number of occurrences of the sequences I am interested in using the seqinr and biostrings packages. I have been using biomaRt to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together.
The issue I am having is that biomaRt returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, Genatlas and NCBI are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?
Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.
Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.
library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences <- function(x) {
  seqMax <- split(x, factor(x$mgi_symbol))
  seqMax <- lapply(seqMax, function(x) max(x$length)) %>%
    unlist() %>%
    as.data.frame()
  seqMax$mgi_symbol <- row.names(seqMax)
  row.names(seqMax) <- NULL
  names(seqMax) <- c("length", "mgi_symbol")
  seqUnique <- merge(seqMax, x)
  seqUnique <- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences <- function(genelist, seqType) {
  ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ <- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ <- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ <- sequ[order(sequ$mgi_symbol), ]
  sequ$length <- sapply(sequ[ ,seqType], nchar)
  sequ <- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs <- function(x) {

  seq <- x %>%
    tolower() %>%
    s2c()

  output <- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXX" = count(seq, 3)["xxx"] %>% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq <- getSequences("Grin1")
targetMotifs <- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable <- data.frame(Length = seq$length)
rownames(targetMotifTable) <- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable <- cbind(targetMotifTable, targetMotifs) %>% as.data.frame()
names(targetMotifTable) <- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb <- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"),
            filters = "mgi_symbol",
            values="Grin1",
            mart = ensembl,
            bmHeader = TRUE)

It would be much easier to help you if you could provide us with a minimal example that reproduces the issue that we can run at home and play with. Especially the specific transcript and the commands you use to get the different numbers of exons.

@terdon I've added my code. Is this a sufficient MWE?

If you have the string you search in X:Y can it be on x+1:y+1? you can look for Sushi, Gviz packages. Also, why don't you simply use grep with your motifs?

@ Llopis Can you clarify what you mean by X and Y? I don't understand. Thank you for the package recs! This might be exactly what I need but it will take me a while to look through it since there are a lot of features available. As for grep... I'm not exactly sure, but I wrote this code last year and I think grep was too slow because I was going over 35000 genes with this code. From what I remember, the count() function had some C wizardry going on that made it faster. I think this is why I used it instead of grep but I can't remember exactly.

@syntonicC For X and Y I mean positions in the genome. The reason I ask this is to know if you are interested in overlapping motifs or not.

@syntonicC One more thing, could you please provide example motifs/sequence and sequences where they should be to play with it? Or point some ids as examples.

@Llopis Sorry, I'm still new to this and can't figure out what you are asking. Does X and Y indicate start and stop positions along a chromosome or something else? Overlapping motifs: If I search for a substring wouldn't this capture any motifs that overlap? I am interested in overlapping motifs. Some motifs I am interested in: "ggcc", "gac", "agga", "gaga". Is this what you are looking for?

@syntonicC Yes, X is the beginning and Y is the end of the region (if this is start or stop of an open reading frame I don't know). Yes, I am interested in this kind of thing, for instance in your data the 3rd and 4th motifs could be overlap in ag*gaga*ta (here in bold and in italics), and some example of gene/transcript (a good model or I will pick one randommly)

@Llopis Ah sorry, I mentioned the gene in my code (and URLs) but did not mention in the question body that it was mouse GRIN1. The coordinates are: 25174720 to 25146695 on the reverse strand. Based on what you said, yes, I do want overlapping motifs as well.

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:
library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)

This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns

I couldn't remember the name of the package ggbio (I was beginning to  doubt it existed). This is far more informative than my answer (+1) BTW, It might be better to overlap the plots of the gene model and the motif, with different colors

I assume you have the sequence in seq1:
# Set the input
seq <- getSequences("GRIN1")
motif <- "ggcc"
# Convert to lower case
seq2 <- tolower(seqinr::s2c(as.character(seq)))
istarts <- seq(from = 1 , to = length(seq2), by = 1)
# Create words of the same size as the motif
oligos <- seq2[istarts]
for (i in 2:nchar(motif)) {
        oligos <- paste(seq2, seq2[istarts + i - 1], sep = "")
}
# Remove exceeding words (words with NA, might need tweaking to motif size
oligos <- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
# Find motif in the sequence
hits <- grep(motif, oligos, ignore.case = TRUE)
#Prepare the plot
positions <- ifelse(1:length(seq2) %in% hits, 1, 0)
plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)

The plot is not the density, but you can choose how to plot it now, as you might be interested in windows sizes of 5 nucleotides or 20.

From here you can iterate through multiple sequences. If you want to find several motifs you can take advantage of the oligos object if they are of the same size.
There is the seqPattern package in Bioconductor for more advanced pattern matching. But surprisingly it doesn't have this simple motif search. I think it would be easier to do this with perl and pyhton than in R.

1: In the code posted the function getSequence returns an error.

I fixed the error in getSequences(). You will need to change your first line of code to seq <- getSequences("Grin1", "cdna") and then also alter this line to seq2 <- tolower(seqinr::s2c(as.character(seq$cdna))), assuming the seqType of interest is "cdna". I think the final line has a typo: the variable is positions not position.

Thank you for your response! However, though I asked for this distribution in my question title, in the question body I stated that I was interested in determining the intron-exon boundaries as well.

Now that I am reading more into this, numbering exons doesn't appear to be that simple. It looks like it is not standardized because there can be internal splice sites within exons as well as exons on different strands. This is probably why I am having so much trouble working this out.

So, is your problem how to calculate the distribution of a motif in a sequence or find an intron-exon boundary? Which answer do you want. Please clarify (I spent quite a bit of time trying to answer this question properly but this format is just for one answer for one question)

Well, it's really both. In the question body I said: "In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms." I should have put this in the question title but it's too late now; I should have been more specific. You answered the question in the title exactly as I asked it so I'm going to accept your answer. Thanks!

Glad to help a bit! It was a nice problem to solve in R.

@syntonicC In case you haven't seen, the other answer is better as it answers both of yours questions.


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

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

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

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

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

The answers can be found on google.


text:
How to best detect the "peaks" in RNA-seq data that are not assigned to any gene?
I encountered that many reads from single-cell RNA seq data were lost in the analysis because not assigned to any gene (genome: galgal6). I am trying to find an approach than could give me all the "peaks" in this data, that is to say regions where there is a high density of reads. Once I get there, I would extract "peaks" that are not assigned to any gene, and get the distance to the closest genes.
I am thinking of using first bedtools genomecov to get the "peaks", and then bedtools intersect to get all reads assigned to a gene. Finally, I extract reads that are not in the output of bedtools intersect to get a new file with all my reads from the peak, that are not assigned to any gene.
Does it make sense this way ? Can you see any other way of doing ?
Here are some tiny data:
> cat reads.bed
chr9    505479  505498
chr9    508014  508037
chr9    514603  514633
chr9    529519  529540
chr9    529519  529540
chr9    529519  529540

> cat tiny.galGal6.chrom.sizes
chr9    24153086

> bedtools genomecov -bg -split -i reads.bed -g tiny.galGal6.chrom.sizes
chr9    505479  505498  1
chr9    508014  508037  1
chr9    514603  514633  1
chr9    529519  529540  3

The last line would correspond to a peak. However, another difficulty is how should I define the thresholds ?
Any help would be more than welcome.

Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use stringTie or similar programs. Ensure that the BAM file you give to stringTie have all of the cells together.

Many pipelines have been set up for this such as software developed for the FANTOM Project at RIKEN. The software for this is available from our website. These were set up for analysis of CAGE sequencing data for exactly this purpose of calling peaks of previously unannotated gene expression. You can also find data for many species (including chicken) from FANTOM5 where this has already been done.


text:
Why PHMM is used For Protein Analysis
I am interested on how profile hidden markov model is used for protein analysis in bioinformatics. Can you guys explain why phmm is more preferable method comparing to simple MSA(Multiple Sequence Alignment) scoring or simple MSA probability scoring? For example we have three MSA sequences as based model and a target to be scored.
MSA Sequences
A B C B C A
A B B C C A
A B C B A A
Target Sequence
A B B B C A
By MSA scoring we can count score for target sequence comparing to a profile MSA. For example above target sequence protein column will counted as match(A), match(B), mismatch(B-majorityMSA(C)), match(B-majorityMSA(B)), match(C-majorityMSA(C)), match(A). Just give score for match and mismatch and we can get total score for a target sequence comparing to a profile MSA sequences.
By MSA probability scoring we can count probability of each protein in each column so for example above
P(target_seq similar to MSA profile) = P(A in MSA sequences)*P(B in MSA sequences)*.... *P(A in MSA sequences) = 1*1*(1/3)*(2/3)*(2/3)*1 = 4/27
Isn't MSA scoring like above examples powerful enough?
Thanks a lot !!

Please [edit] your question and tell us exactly what type of analysis you are thinking of. By MSA, I assume you mean multiple sequence alignment (or are you thinking of some sort of mass spectrometry?), but I don't see how the two are related. Are you thinking of HMMs of protein domains? Something else? Are you referring to finding protein homologs? Tools like PSI-BLAST or HMMER?

I've edited my question, yes it's all about protein sequence analysis.

The simple probabilistic scoring scheme you describe to query a target sequence against a (protein) family represented as a multiple sequence alignment (MSA) is actually similar in essence to a Position Weight Matrix (PWM), although PWMs are more refined. PWMs are powerful rapid tools often used for motif analyses.
However, your approach (and PWMs) don't handle insertions or deletions (indels), and you can't easily add these into your model (for example, consider that 2 1bp indels should not score the same as 1 2bp indel, and indels could occur within your MSA as well as in your target sequence). You could try and add some arbitrary gap opening and extension penalties somehow (like BLAST) but then you'd abandon your probabilistic model! By contrast, pHMMs have a fully fledged probabilistic indel model with insertion and deletion states.
Finally, of practical note, there are some excellent implementations of pHMMs for protein sequence analysis, such as HMMER3 and HHSuite.

I am sorry, still don't get it. What do you mean by 2 1bp indel and 1 2bp indel ? Can you give some illustration?

@RamandikaPranamulia 1 2bp deletion: ATG--AT and 2 1bp deletion AT-GA-T

I recommend you to read this paper in which the applications of HMM in protein sequence alignments are discussed.