text:
How to validate that BAMs have been downloaded correctly?
I currently have several hundred BAM files which were downloaded by someone else. These have remained untouched---before working with them, I would like to double-check that these BAMs have been fully downloaded.
I don't MD5 Checksums to look at.
In order to validate the BAMs, I would use something like Picard's ValidateSamFile, https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile
java -jar picard.jar ValidateSamFile \
      I=input.bam \
      MODE=SUMMARY

I could write a bash script, and run through these BAMs, and then read through the Picard output. That is checking whether the BAM is valid though, not whether it was downloaded fully (e.g. the BAM could be downloaded correctly, but be inherently flawed in some other way).
If the BAM is downloaded fully, isn't there some "end of file" bytes I could quickly check? Are there other more efficient methods?

samtools quickcheck is all you need. From the manual:

Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only).
Data in the middle of the file is not read since that would be much more time consuming, so please note that this command will not detect internal corruption, but is useful for testing that files are not truncated before performing more intensive tasks on them.
This command will exit with a non-zero exit code if any input files don't have a valid header or are missing an EOF block. Otherwise it will exit successfully (with a zero exit code).

Just a note if your cluster has older modules, or old binaries in your search path: samtools quickcheck requires samtools ≥ 1.3. For older versions you can samtools stats BAMFILE |& grep -q 'EOF marker is absent' && whatever, but this still tries to compute stats for the non-corrupted parts of the BAM so it's not the fastest solution.


text:
Output of Seurat FindAllMarkers parameters
I compared two manually defined clusters using Seurat package function FindAllMarkers and got the output:

Now, I am confused about three things:

What are pct.1 and pct.2?
How come p-adjusted values equal to 1? What does it mean?
If we take first row, what does avg_logFC value of -1.35264 mean when we have cluster 0 in the cluster column? Is it that in cluster 0 the Cttnbp2 gene downregulated by a factor of 2^1.35264?

https://satijalab.org/seurat/de_vignette.html

You have a few questions (like this one) that could have been answered with some simple googling. I suggest you try that first before posting here.

I could not find it, that's why I posted. Thank you @heathobrien! Should I remove the Q?

pct.1– The percentage of cells where the gene is detected in the first group

p_val_adj– Adjusted p-value, based on bonferroni correction using all genes in the dataset.

This is not also known as a false discovery rate (FDR) adjusted p-value.

An adjusted p-value of 1.00 means that after correcting for multiple testing, there is a 100% chance that the result (the logFC here) is due to chance.

You are correct. The expression of Cttnbp2 donwn-regulated by that factor. If you were to look in the table at the cluster 1 DE genes you would see Cttnbp2 with a logFC value of 1.35264. When you are looking at genes differentailly expressed in only two clusters, the sign is the only thing that changes.

Vignette posted above by @heathobrien.

How the adjusted p-value is computed depends on on the method used (use.test). E.g. "DESeq2" uses Independent Hypothesis Weighting. While Bonferroni correction does indeed control the FDR as a side effect it is build to control the family-wise error error rate. An adjusted p-value of 1 does not mean that the result is due to chance, it just means that it cannot be excluded (by the methods used) that such a difference would be observed was there no true difference in expression.

source: Stack Exchange

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

sort: votes


text:
Expression of a gene in different groups
I would like to check the expression of a gene in different groups like Disease vs Normal samples. I want to make a plot out of that to check whether it is significant or not.
From this paper lncRNA I see that Figure 1C they used RPKM value. But, I'm using read counts data. My question is can I use log cpm to make a plot like that or should I need to use only RPKM?

To check if a gene is significantly differential expressed between groups, is not done with a plot but with linear models (e.g., limma trend or voom, edgeR, or DEseq2). Doing a t-test with RPKM data like in the paper is certainly not recommended (but reviewers are usually no statisticians or bioinformaticians, so that unsound statistics is still published). My advice is to keep count data for statistics in limma (or other tools I mention), and make violin plots with cpm.

log-scale is very important for differential expression analysis but I think why we do this is best addressed as a separate question.

You should never use RPKM. It’s simply obsolete in the age of paired-end sequencing, and has been replaced by FPKM (which is, strictly speaking, a synonym).
The linked blog post explains more generally the problems that measures such as FPKM and CPM suffer from. A more robust measure is the TPM (transcripts per million), which scales CPM by the lengths of the transcripts.
Lastly, none of these measures perform very well for comparing gene expression across replicates. Since you have read count data, you should run DESeq2 or edgeR to perform differential expression analysis across all genes, and go from there.
As the comment noted, the paper you cite does not represent scientific best practice. In fact, a competent reviewer would have rejected this analysis.

For visualization purposes, using log(cpm) is fine. But plot don't check if differences are significant or not, statistical tests do. You can certainly add the results of a proper statistical test to a plot - like adding asterisks to denote a significant difference, for example.

You can use the raw counts to normalise your data. It is import to do so before doing a differential expression analysis: otherwise you will detect more reads from samples with higher coverage and from longer genes. A gene with longer mRNA molecules will produce more reads than a shorter one for the same number of transcript molecules. Samples will not produce the same number of reads due to variation in cDNA extraction, amplification, library preparation, etc.
This does not inherently mean that a sample expresses every gene more than another sample. Hence we compare samples by relative rather than absolute expression (for microarray and RNA-Seq data). This is in principle, very similar to using housekeeper genes in qPCR experiments except we adjust for total RNA or reads, gene length, and batch effects.
There are many ways to achieve this but I recommend using the limma package in R. This provides the “voom” normalisation procedure. “voom” is documented in a separate follow-up paper to the limma paper. It produces normalised data (on a log-scale) which is compatible with downstream analysis provided by the package (including built-in Venn diagram and volcano plot functions). See the package vignette for more details.
Of course, everything here only applies to bulk RNASeq experiments and procedures for single-cell analysis are different.

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


text:
Demultiplexing and preprocessing for custom single-nucleus Drop-seq data
I am trying to reproduce the preprocessing of paired-end sequencing reads in Lake et al. (ref. 1).

First, paired-end reads were filtered out if read 1 had more than four non-T bases in the last ten bases (to remove all non-poly(T)-captured contaminated reads), or had one or more bases with a poor quality score (<10). Cell barcode and UMI information were then inferred from the first 12 bases and the next 8 bases of read 1, respectively. The right mate of each read pair was trimmed to remove any portion of the SMART adaptor sequence or large stretches of poly(A) tails (6 consecutive bp or larger).

I was able to filter out reads as described in the first sentence.
However, I am unsure on which tool or procedure to use to demultiplex based on cell barcode and UMI, and to trim the right mate. Do you know what I can use or something I can read to understand how to do it?

Lake, B. B. et al. Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nature Biotechnology 36, 70–80 (2017).

I didn't notice the sentence:

Paired-end sequencing reads were processed largely as described (http://mccarrolllab.com/wp-content/uploads/2016/03/Drop-seqAlignmentCookbookv1.2Jan2016.pdf ), with additional correction steps.

This should provide instructions on how to lead the preprocessing and mapping of data.


text:
What is cellranger doing in comparison to other methods?
I've recently started working with the 10X-Genomics platform with Illumina (MiSeq and HiSeq) for single-cell RNA-Seq. I've been recommended the "cellranger" (version 2.1.0) which I understand handles the barcoding of the platform and performs processing of the data. For example:
1) convert raw Illumina .bcl files into .fastq (demultiplex)
cellranger mkfastq

2) perform alignment, barcode counting, and generate gene-barcode-martices
cellranger count

However, it is unclear which underlying bioinformatics procedure (with which parameters) are being used here. For this pipeline, I wish to know how the trimming, alignment and counts are being performed. This is important to understand whether further normalisation and quality checking is necessary. What steps would I need to do to perform a similar analysis?
If for example, I wish to document in more detail what has been performed here in a reproducible manner or reproduce the analysis pipeline on a different platform (such as Nadia) to compare them.
Ideally, I would like to know the analogous pipeline using open-source tools but at the moment I am first seeking to understand what is going on under the hood when I am running the cellranger.

Have you looked at the cellranger source code (e.g., the alignment command is generated here)? If you want the exact commands run that's the easiest method.

Thanks for the link to the source code. That is helpful but full for the benefit of those not well-versed in Python (or these tools) more explanation of what these parameters are doing would be helpful.

As an update: kallisto | bustools, dropEst, and DropSeqPipe are valid alternatives for this task.

cellranger mkfastq is not necessary anymore.  It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself.  Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them.
Cellranger count aligns the reads, filters away duplicates based on UMIs, tries to figure out which cell barcodes really captured cells, and tells you for each cell barcode how many reads hit each gene.  You can go through the .bams, and see which reads were counted for what gene and what barcode, based on the flags and other tags in the .bam
If you want another pipeline to compare to, you can use the McCarroll lab's DropSeq pipeline, which is a bunch of little scripts, around the output of STAR aligner (just like 10XGenomics uses).  I have put the same data sets through both 10X and Dropseq, and they return about the same thing.

Great answer but it would help to have more detail on how the cellranger tool is performing the alignment, is this with the default parameters for STAR and samtools? I've been advised to "demultiplex" MiSeq output but that may be for older versions of cellranger.

The best approach would be to contact 10x Genomics support: support@10xgenomics.com . They usually respond within a few hours.
Otherwise, you can check the methods from the first 10x Genomics paper:

The Cell Ranger Single-Cell Software Suite was used to perform sample
  demultiplexing, barcode processing and single-cell 3′ gene counting
  (http://software.10xgenomics.com/single-cell/overview/welcome). First,
  sample demultiplexing was performed based on the 8 bp sample index
  read to generate FASTQs for the Read1 and Read2 paired-end reads, as
  well as the 14 bp GemCode barcode. Ten basepair UMI tags were
  extracted from Read2 (14 libraries were made with 5 bp UMI tags, as
  noted in Supplementary Table 1, due to an earlier iteration of the
  methods. For these samples, 5 bp UMI tags were extracted from Read2.).
  Then, Read1, which contains the cDNA insert, was aligned to an
  appropriate reference genome using STAR35. For mouse cells, mm10 was
  used and for human cells, hg19 was used. For samples with mouse and
  human cell mixtures, the union of hg19 and mm10 were used. For ERCC
  samples, ERCC reference
  (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt)
  was used.
Next, GemCode barcodes and UMIs were filtered. All of the known listed
  of barcodes that are 1-Hamming-distance away from an observed barcode
  are considered. Then, the posterior probability that the observed
  barcode was produced by a sequencing error is computed, given the base
  qualities of the observed barcode and the prior probability of
  observing the candidate barcode (taken from the overall barcode count
  distribution). If the posterior probability for any candidate barcode
  is at least 0.975, then the barcode is corrected to the candidate
  barcode with the highest posterior probability. If all candidate
  sequences are equally probable, then the one appearing first by
  lexical order is picked.
UMIs with sequencing quality score >10 were considered valid if they
  were not homopolymers. Qual=10 implies 90% base call accuracy. A UMI
  that is 1-Hamming-distance away from another UMI (with more reads) for
  the same cell barcode and gene is corrected to the UMI with more
  reads. This approach is nearly identical to that in Jaitin et al.4,
  and is similar to that in Klein et al.8 (although Klein et al.8 also
  used UMIs to resolve multimapped reads, which was not implemented
  here).
Last, PCR duplicates were marked if two sets of read pairs shared a
  barcode sequence, a UMI tag, and a gene ID (Ensembl GTFs GRCh37.82,
  ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz
  and GRCm38.84,
  ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz,
  were used). Only confidently mapped (MAPQ=255), non-PCR duplicates
  with valid barcodes and UMIs were used to generate gene-barcode
  matrix.
Cell barcodes were determined based on distribution of UMI counts. All
  top barcodes within the same order of magnitude (>10% of the top nth
  barcode, where n is 1% of the expected recovered cell count) were
  considered cell barcodes. Number of reads that provide meaningful
  information is calculated as the product of four metrics: (1) valid
  barcodes; (2) valid UMI; (3) associated with a cell barcode; and (4)
  confidently mapped to exons.

Of course, they have since updated the library prep and the Cell Ranger software, so the exact protocol is somewhat different.


text:
update dbSNP ID
I have used dbSNP build 138 to tag 'rs' ID in my VCF file. But later I came to know that now dbSNP build 150 has been released. I want to update the rs IDs in my VCF file.
I tried but I found a problem for which I need suggestion. With dbSNP build 138, I am getting 450076 total variants, 248236 as known variants and the rest as novel. But when I use the latest dbSNP file, All_20170710.vcf.gz, I am getting same total variants but 0 as known and all as novel which is definitely wrong.
When I checked the reason for this, I found that in the All_20170710.vcf.gz VCF file, all chromosomes are written as numbers like 1, 2, instead of chr1, chr2 while in my vcf file, all chromosomes are written as chr1, chr2. The dbSNP build 138 file also has the chromosome names with the chr prefix, and that's why 138th build was giving good results.
One solution would be add the chr word before chromosome number. Is that an efficient way? How do I do that? If not, what is the correct way to use  the latest dbSNP All_20170710.vcf.gz to tag rs IDs in a vcf file which contains ~32 Billion SNPs instead of dbsnp_138.vcf which contains only ~6 Billion SNPs? Any suggestions??
Thanks

It sounds like you downloaded the dbSNP release for the GRCh38 or GRCh37 genome releases from the Genome Reference Consortium (GRC) which name the chromosomes without the chr prefix and you were previously using a file that had been built using either hg19 or hg38 which use the chr.
As far as I know, dbSNP only offers vcf files built against the GRC genome builds, so you (or someone else) had probably already modified the old file you were using and changed the chromosome names.
Also note that if you want the latest version, you should download ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/All_20180418.vcf.gz.
Anyway, yes, you can simply change the VCF to use chrN instead of N. Note, however, that you will also need to change the name of the mitochondrial chromosome from MT to chrM
You can do all this very easily with sed:
zcat All_20170710.vcf.gz |
    sed -E '/^[^#]/{s/^([^#]+)/chr\1/; s/^chrMT/chrM/}' > All_20170710.fixed.vcf

The command above will work on most sed implementations, including Linux (GNU sed) and also macOS (BSD sed).

@heathobrien, @terdon : A little clarification here. I have downloded dbsnp_138.hg19.vcf.gz which use chr prefix and in file All_20170710 from GRC contains number only.

@Lot_to_learn whoops, yes sorry. GRC uses numbers only, it's the hg19 and hg38 that use chr. Answer edited. Also, where are you downloading these files from?

@Lot_to_learn please see the updated answer. The first version would not have corrected the names for chrX and chrY.

I have found more rs IDs from new dbSNP as comparison to old one. For one patient I got 3200 approx new rs IDs. Although I was expecting more new rs IDs. Thanks. This is working well.

@Lot_to_learn did you use the very latest edited version of my answer? The first had a couple of problems. To check, make sure that you have chrMT, chrY and chrX in the file.

A little correction in my previous comment. I don't know why, when I checked no. of rs ID's using grep, it shows big number but gatk tool shows very less known rs ID.So when I use latest version of dbSNP I found very less known and most of the mutations are novel according to gatk variantEval. When I checked, I found that its reference is GRC38. So I use 2017 version of dbSNP. Now I am getting lot more known variants in variantEval. Now I am OK with 2017th version of dbSNP. Thanks.

Keep in mind that the different names sometimes correspond to different genome builds, so you should make absolutely sure that you are using the right (compatible) reference genome build before doing this. Eg you can't use a dbsnp set from b37 with a callset made using hg38 and vice versa.


text:
How can I easily get the read size distribution of reads mapping on a certain set of regions?
Suppose I have a BAM file indicating where reads in a library have mapped, and a bed file describing a set of genomic regions.
Is there a way to easily get the size distribution of the reads mapping on this set of regions?

This can be done rather simply by combining this samtools based answer and bioawk:
samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1

The samtools command extracts the alignements in the region given with the -L option, in sam format. The -c sam option of bioawk makes the $seq variable (corresponding to the read sequence field) available in the awk command.
The final sort ensures the size classes are ordered numerically.
We obtain a two-column size    count table.
Examples with timings:
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s

Since sam format is actually simple to parse, one can try to use mawk instead for slightly higher speed (the sequence being in the 10th field):
samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1

(Contrary to bioawk, the default output field separator is not a tabulation.)
Example:
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s

Via BEDOPS bedmap and bam2bed, along with the hist functionality in bashplotlib
$ bedmap --echo-map-size --multidelim '\n' intervals.bed <(bam2bed < reads.bam) | hist

 77|           o
 73|           o
 69|           o
 65|           o
 61|           o
 57|           o
 53|           o
 49|           o
 45|           o
 41|           o
 37|           o
 33|           o
 29|           o
 25|           o
 21|           o
 17|          oo
 13|         ooo
  9|         ooo
  5|         ooo
  1| o    o oooo
    -----------

------------------------------
|          Summary           |
------------------------------
|     observations: 118      |
|    min value: 16.000000    |
|      mean : 34.415254      |
|    max value: 36.000000    |
------------------------------

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


text:
Convert R RNA-seq data object to a Python object
I have done some work in R and would like to try a Python tool.
What is a good way to import the data (and its annotations etc) as a Python object?
I am particularly interested in converting a Seurat object into an AnnData object. (Either directly or as a Python object which can be converted into an AnnData.)

A simple solution for converting Seurat objects into AnnData, as described in this vignette:
library(reticulate)
seuratobject_ad <- Convert(from=seuratobject, to="anndata", filename="seuratobject.h5ad")

Alternatively, one can use Loom, "a file format designed to store and work with single-cell RNA-seq data".
In R, install LoomR:
devtools::install_github(repo = "mojaveazure/loomR")

Convert from Seurat object to loom:
pfile <- Convert(from = pbmc_small, to = "loom", filename = "pbmc_small.loom")
pfile$close() # close loom objects when finished using them.

Then import loom object in Python using loompy, or directly as AnnData:
scanpy.api.read_loom

Alternatively, see feather.

Or export as some text format (csv, json) then import into Python.

I would recommend to follow the KISS rule. If you can export the data in .csv or json better, but that depends on what is the original data format

The accepted solution is probably the best for older objects of type seurat created with Seurat package v2. For newer Seurat Objects, there is a new tool designed specifically for this purpose, called SeuratDisk.
There is a nicely documented vignette about the Seurat <-> AnnData conversion. In short:

In R, save the Seurat object as an h5Seurat file:

library(Seurat)
library(SeuratDisk)

SaveH5Seurat(SeuratObject, filename = "file.h5Seurat")
Convert("file.h5Seurat", dest = "h5ad")

In python, simply load the file:

import scanpy

adata = scanpy.read_h5ad("file.h5Seurat")

Since both Seurat and Scanpy are evolving rapidly, the safest way is to
save your Seurat object into MatrixMarket format + metadata csv, and
import them into Scanpy.

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


text:
How to download gene expression data from NCBI gene database
In the NCBI gene database, I can add the expression tracks (circled in picture blow) through 'Tracks' button, but How I can download the expression data directly, not just look the picture?

You can download raw data from SRA (NCBI). On this website under RNA-Seq alignments, you'll find the samples. Click on them, for example SAMN05231885, and then click on PRJNA325427 to see the link to the SRA data (click on 1). Here you can download the raw fastq files of the experiment.


text:
Subset data frame by similar elements from rows and column between two data frames
My first data frame whose first column i intend to use which is PATIENT_ID
PATIENT_ID  SEX RACE    FAB AGE BM_BLAST_PERCENTAGE WBC
TCGA-AB-2805    Male    WHITE   M0  77  67  92
TCGA-AB-2822    Male    WHITE   M0  65  99  2.9
TCGA-AB-2831    Male    WHITE   M0  59  61  1.4
TCGA-AB-2890    Male    NH/C    M0  69  60  6.1
TCGA-AB-2907    Male    WHITE   M0  69  82  1.6
TCGA-AB-2913    Male    WHITE   M0  61  95  61.6
TCGA-AB-2922    Male    WHITE   M0  83  75  77.9
TCGA-AB-2936    Female  NH/C    M0  60  50  2.1

My second data frame where i have the respective patient ID as individual column with expression values
Hugo_Symbol TCGA-AB-2803    TCGA-AB-2805    TCGA-AB-2806    TCGA-AB-2807    TCGA-AB-2808
A1BG    0   0   0   0   0
A1CF    0   0   0   0   0
A2BP1   71.4095 48.5475 50.8339 45.7498 47.1328
A2LD1   2.0576  3.3467  5.7288  0.993   5.4988
A2ML1   467.0782    16.0643 7.0019  343.5948    255.9152
A2M 11.3169 0.6693  17.823  116.1867    182.2467
A4GALT  0   2.6774  1.2731  4.9652  0.7855
A4GNT   0   0   0   0   0
AAA1    514.4033    532.7979    504.1375    697.1202    443.8335
AAAS    0   0   0   0   0
AACSL   190.3292    74.2972 84.6595 559.0864    317.3606
AACS    0   0   0   1.9861  0
AADACL2 0   0   0   0   0
AADACL3 0   0   0   0   0
AADACL4 1.0288  0   0.6365  0   0

What I tried so far transposing the second data frame so that i can get the individual column as my first column ,and subset by using dplyr library but so far im not able to ..any help or suggestion how to do it would be highly appreciable
I would like to have my expected data frame as such
PATIENT_ID  Gene B  GeneA   GeneC   GeneA   GeneA
TCGA-AB-2803-03 0   0   0   92.5926 876.5432
TCGA-AB-2805-03 0.6693  1.3387  0   92.6305 421.6867
TCGA-AB-2806-03 0   14.1884 16.3654 119.599 625.0796
TCGA-AB-2807-03 2.006   9.3247  13.5154 178.4211    537.2393
TCGA-AB-2808-03 0   42.3016 28.3975 170.2042    653.5742
TCGA-AB-2810-03 0   1.5668  8.9665  152.0747    1087.5648
TCGA-AB-2811-03 0   0   6.9556  30.3103 1414.6602
TCGA-AB-2812-03 0   14.3289 11.9466 120.3591    962.5575
TCGA-AB-2813-03 0   3.4938  1.3813  51.1883 594.7593
TCGA-AB-2814-03 0   56.8659 118.7051    87.1817 472.6912
TCGA-AB-2815-03 0   2.0692  3.6986  104.6575    502.5234

okay so I what im expecting is I would like to have total 18 patient which are in first data frame my patient ID so I would like to use the patient ID and  use the same to subset the 18 patient from my second data frame with their expression values my second data frame 22000174 in dimension .so my expected data fram would be 2200018 , so i need my rows as my genes and columns as my patient sample ID..

Are you sure that those numbers on that expression matrix are hugo symbols?? I highly doubt it. BTW I think this is not about bioinformatics...

well i just put it as since it didn't contain any symbol but entrez ID

If you're still looking for an answer on this, I would use the gather and spread functions in the tidyr package (part of the tidyverse suite of packages - which also includes dplyr and ggplot2:
library(tidyverse)
df2_long <- gather(PATIENT_ID, value, -Hugo_Symbol)
df1 %>%
  left_join(df2_long) %>%   # merges the 2 data frames together
  spread(Hugo_Symbol, value)

You can take the PATIENT_ID from the first dataframe and use them to select from your second dataframe. Than transpose, so something like this:
ID_Names <- First_df$PATIENT_ID

rownames(Second_df) <- Second_df$Hugo_Symbol

Selected_df <- Second_df[, ID_Names]

Expected_df <- t(Selected_df)

well i m running into an error "Error: Columns TCGA-AB-2831, TCGA-AB-2907, TCGA-AB-2922 not found" may be one problem i think is my ID_Names vector is in character , but the column i m trying to map are in numeric ...any suggestion

Are you sure these three patient_IDs are in the second file? You select from the header which is not numeric.

i mean all my header in the second data frame are numeric...so do i need to convert them to character if yes how do i do that i tried converting them to character it changes my data frame

Can you show: class(colnames(second_df)), I see characters in your colnames, so I doubt that you are right about numeric header. Are you sure every ID is in the colnames of your second data frame?

The error is pretty clear, you miss the three IDs in your second df. Check the colnames and look for the specific three IDs in the error.

[1] "character" my first colname is character as i see and rest as numeric when i hover my cursor over the rest

yeah now its solved..,i think those no ID in my second df, but now i transposed the data frame but the gene names are not coming..

Make rownames with the gene names for the second df, before you select. Something like this: rownames(second_df) <- second_df$Hugo_Symbol. If you then select the columns the gene names are stored as rownames.


text:
How to get results from Homo.sapiens package in bioconductor for a specific reference
I want to use the Homo.sapiens package in Bioconductor to retrieve the chromosome location start and end for each gene symbol in a specific reference (e.g. hg19 or hg38).
Right now, I am using the following code:
> select(Homo.sapiens, keys=c('KRAS'), cols <-c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype='SYMBOL')

  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863

My question is two-fold:

How can I specify the human genome reference I want to retrieve the data from? (e.g. hg19 or hg38)
Why do I receive multiple and different mappings through Homo.sapiens for some genes? If I check the locations through GeneCards I get:

chr12:25,204,789-25,250,936 (GRCh38/hg38)
chr12:25,357,723-25,403,870 (GRCh37/hg19)

I couldn't find how is this package build, my guess is that is as in the example of the OrganismDBI package. See that vignette section to build your own data package with the genome you want.
The locations are different because you ask for the transcript start and end, if you add the "TXNAME", you will (probably) see there are different transcripts, so different mappings are reported (untested).
Update
The example of the OrganismDBI package builds correctly the Homo.sapiens package, which appears to be linked with hg19. I provide here the relevant code to build the package for hg19 and hg38.
> source("https://bioconductor.org/biocLite.R")
> biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
> biocLite("TxDb.Hsapiens.UCSC.hg38.knownGene")

> library(OrganismDbi)

# instructions for hg19, for hg38 change accordingly
> gd <- list(join1 = c(GO.db="GOID", org.Hs.eg.db="GO"),
             join2 = c(org.Hs.eg.db = "ENTREZID",
                       TxDb.Hsapiens.UCSC.hg19.knownGene = "GENEID"))
> destination <- tempfile()
> dir.create(destination)
> makeOrganismPackage(pkgname = "Homo.sapiens.hg19", graphData = gd,
                      organism = "Homo sapiens", version = "1.0.0",
                      maintainer = "Maintainer<maintainer@email>",
                      author = "Author Name", destDir = destination,
                      license = "Artistic-2.0")

Creating package in /var/folders/sf/tfrswpc96jl8b4520mn784440000gn/T//RtmpzR4YnY/file111b490af70f/Homo.sapiens.hg19

# installing package
> install.packages(destination, repos = NULL, type="source")

Then, you can load and use the packages with:
> library(Homo.sapiens.hg19)
> library(Homo.sapiens.hg38)

> select(Homo.sapiens.hg19, keys=c('KRAS'),
     cols = c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'),
     keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863
> select(Homo.sapiens.hg38, keys = c('KRAS'), cols =
       c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25204789 25250931
2   KRAS   chr12 25209168 25250936
3   KRAS   chr12 25209431 25250803
4   KRAS   chr12 25233819 25250929

You're right about TXNAME, I'll check about building my own data package, even if I'd like a tested implementation. Thanks

@gc5 that is the tested implementation, this data packages only collect the information available from other packages. So the data isn't altered

@gc5 Thanks for adding the code you used. I only tweaked the format

Yes, thanks. I still don't understand why the hg19 package does not provide a coherent mapping of transcript start and end with GeneCards, but I think that's minor.


text:
Books on bioinformatics algorithms
I'm looking for a book about bioinformatics algorithms, such as alignment, BLAST search, and variant calling.
I'm hoping reading about this subject will give me a deeper understanding of the foundations of bioinformatics, and I'm also interested in developing my own tools and algorithms for sequence and variant analysis. I studied applied mathematics and computer science at university, so I don't shy away from math-heavy books.
I've listed some candidates below.

Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleid Acids, 1st edition, Durbin et al. 1998
An Introduction to Bioinformatics Algorithms, 1st edition, Jones et al. 2004
Algorithms in Bioinformatics: A Practical Introduction, Sung 2010

Durbin et al. is quite old, Sung is the newest.
I am looking for books with:

material (algorithms) on alignment, assembly, and variant calling if possible.
New books are preferred

"Biological Sequence Analysis" is excellent, you say it's old but it's been re-printed many times with edits and is still very relevant. My edition is the 10 re-reprinting from 2010, they may be an even more recent one

Ziheng Yang's "Computational Molecular Evolution" is also good I recall

@Llopis it's common on biology.stackexchange.com to ask for book recommendations. I asked the same question there, and was recommended that I ask it here instead. For what it's worth, I think the question is better suited here, as it is more likely to reach relevant people here.

The asker is looking for books that meet various fact-based criteria — it covers topics A, B, C, etc. A book recommendation objectively either meets or does not meet those criteria. If it does not meet the criteria, that is a problem with the answer, not the question. It's difficult to see what is open-ended or opinion-based about meeting plainly-stated criteria. Voting to reopen.

Completely agree with Alex, voted to reopen.

Here are a couple books I'd recommend:
Dan Gusfield's Algorithms on Strings, Trees, and Sequences is a deep and wide treatment of aligning, searching, and processing strings, trees, and sequences.
Warren Ewens' Statistical Methods in Bioinformatics devotes a chapter to BLAST and the math underneath it.
Edit - Another book that may be useful for some mathematical underpinnings would be Michael S. Waterman's Introduction to Computational Biology.

I agree that algorithm books, like Gusfield's, are hugely important. Personally though, I have that covered decently already.

The Waterman book looks quite interesting, it has some material that I haven't seen in other books.

If you're interested in variant calling, I just co-authored an O'Reilly book called Genomics in the Cloud that covers the GATK Best Practices from a scientific standpoint (germline and somatic short variants + somatic copy number alterations), as well as technical considerations like how to run it efficiently at scale (that's where the cloud part comes in).
You can read more about it on the GATK blog; there's an announcement and a second more personal blog post that I wrote to explain how I approached the topic. Happy to discuss more if you're interested.
Oh and you can preview the first few chapters for free on Amazon if you like.

Bioinformatics and Functional Genomics by Jonathan Pevsner is a good one, he goes over the biology as well as how the algorithms work, and provides real-world examples.

This one is very nice, a colleague of mine happens to have it. It doesn't seem to go into much detail with the math and algorithms though (which is not a deal breaker or anything), probably because it's targeted mostly at biology students. I think it has a very nice way of tying the algorithms to the biological context. Perhaps this book is a good introduction, and to get more mathematical details one can read papers on BLAST, aligners such as BWA, and so on.


text:
A good tool for gene locus visualization
Does anyone know of a good tool for visualizing gene loci? Basically, I'd like to show how a group of genes look in different strains of a bacteria. I have annotation of course, I just need a tool to plug it into and get an image of the locus to compare.
I'm talking something basic like this:

I ended up using Easyfig, which is made for the task and requires little tinkering around, once fed the genbank files.
(Granted, being so simple, it's hard to make the image more advanced, but it serves my purposes for now).

Here's an example of my image.

Could you post an image from the software (to help future readers) ? Thanks for coming back

I think in R, Gviz is the best choice. I have used it for human sequences but it will probably also be suitable for Prokaryotes.
Another option is MS powerpoint, if you only need the arrows...

yah +1 gviz is great, i've used it for bacteria

I wanted to avoid putting something together manually with powerpoint, even though I suspect that's how my supervisor does it...
I'll take a look at Gviz though!

GFF2PS was developed for exactly this sort of thing. It takes a GFF file as input, and will produce a pretty PS output showing the annotations. It was used to create the poster for Celera's human genome draft [1].
Here's a very simple example of what it can do:

And a more complex one from the mosquito genome paper [2] (click here for a full size, 1.5M gzipped PostScript file):

If you are happy with the defaults, it is very easy to use, but otherwise it is a bit of a pain to configure. It is absurdly powerful though.

Venter et al. Science 291(5507):1304-1351 (2001)
Holt et al. Science 298(5591):129-149 (2002)

gggenes is a ggplot2 based solution and looks very appealing. Plus, it has a nice introduction with examples.
library(ggplot2)
library(gggenes)
ggplot(example_genes, aes(xmin = start, xmax = end, y = molecule, fill = gene)) +
  geom_gene_arrow() +
  facet_wrap(~ molecule, scales = "free", ncol = 1) +
  scale_fill_brewer(palette = "Set3") +
  theme_genes()


text:
state-of-the-art methods to obtain SNVs from circulating tumor DNA?
What is the state-of-the-art to obtain SNVs from circulating tumor DNA after Illumina sequencing? This is, finding a list of nucleotide locations that show tumour SNVs filtered out of any germline polymorphisms.
Found so far:

SiNVICT: ultra-sensitive detection of single nucleotide variants and
indels in circulating tumour DNA.
OpenGene/ctdna-pipeline
MAGERI: article
MuTect: running it against a synthetic normal for unmatched cfDNA samples that don't have normal, then filtering for 1000G germline mutations.
CaVEMan: running it against a reference normal for unmatched cfDNA samples that don't have equivalent normal, then filtering for 1000G germline mutations.
ShearWater: corrects against germline mutations by training against a set of healthy normal samples.

Do you think you can narrow a bit more? In what aspects are you interested? What tool(s) have you tried?

@Llopis I have now extended the question and added 3 more pieces of software that I found relevant.

Yes, precisely because you can grow the list ad infinitum I wanted you to narrow the question. In general those open ended questions don't work well in the StackExchange format because they don't have an answerable with a correct answer. And tend to results in a list of software, that although interesting and useful it might fit better for a blog or other type of resources. That's why I asked if you have some problem with any of these tools.

The problem is that, presumably, all of these are "state of the art", so choosing one will be largely a question of opinion. Therefore, this is a very interesting discussion, but not really something that can be given a definitive answer of the sort that the Stack Exchange model requires.


text:
GOseq analysis with evidence code filter
When using GOseq analysis on RNA-seq data, I often find many 'false positives'. What I mean with that is that some genes in certain categories are not really involved in the process, but are only involved due to 'Electronic Annotation' (IEA). Most of the genes are IEA. A list of all evidence codes is here.
Is it possible to (easily) do GOseq analysis without this evidence group?
I mean, I know how to do it the hard way (make a custom GO-genes annotation in e.g., biomaRt removing all IEA evidence code genes, and perform GOseq with custom genome). But is there also a easy way in GOseq available?

For what it is worth, there have been papers that say that the quality of the IEA is not so bad (as bad as other categories), also the GO says that the evidence codes are not designed to be used as quality filters.

Thanks @Llopis, I mostly agree with their statement about IEA quality. However, I was thinking of ways to prevent getting enriched GO terms such as 'spermatogenesis' in female mice experiments... Another motivation is to get 'real' target genes for further investigation (often you get false positives when further investigating).

Then I recommend to give a look at the topGO package, it should filter most of them out, via taking into account the DAG structure of GO.

@Llopis, do you mean with elim function? That reduces the number of significant GO terms, but won't eliminate false positive genes (targets) for further research.

Yes, I mean the elim and the weight01 methods. I agree that they won't eliminate false positive genes but it will let you focus on genes  on categories that are really relevant. Anyway, let's see if someone came up with a solution, even if it is the long one.

I think the answer to the quesiton is that there is not easy way to do it. But as a comment i'd say that I would never trust a GO analysis to provide "answers", only suggestions for further analysis, and as such any genes in significant categories are always going to have to be examined by a domain expert.

@Llopis, can you direct me to the papers about the quality of IEA and about not using evidence codes for filtering. I am interested in their motivation, thanks.

@b.nota sorry I though I had the references but I don't. I am looking for those articles, if I find them I'll post in a new question and let you know

@b.nota I found the reference! I found that Gene Ontology semantic similarity tools: survey on features and challenges for biological knowledge discovery recommends (under the section "Influence of GO evidence codes") using all the data and evidence codes, because the tools (that make these associations IEA) are getting better to assign a term (and it provides references to semantic similarities with IEA and without IEA). Edit: Now perhaps I should write a question on itself...

@Llopis, thanks for getting back with the reference.

A bit weird to answer my own question, but for who's interested I can show the 'hard' way using biomaRt and then use goseq with gene2cat defined.
So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.
library(biomaRt)

table <- read.table(file, header = T, stringsAsFactors = F)

ensembl <- useMart("ENSEMBL_MART_ENSEMBL")

ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)

allGenes <- table$Geneid

GOs <- getBM(attributes=c('ensembl_gene_id', 'go_id', 'name_1006', 'go_linkage_type', 'namespace_1003'),
      filters = 'ensembl_gene_id',
      values = allGenes,
      mart = ensembl)

dim(GOs)
[1] 477082      5

table(GOs$go_linkage_type)
          EXP    IBA     IC    IDA    IEA    IEP    IGI    IMP    IPI    ISA    ISM    ISO    ISS    NAS     ND    RCA    TAS
 36931     36  27918    525  34038 210114    630   6891  30219   6205   4234     20  77365  24673    609  12821     99   3754

# I only want BP and no IEA evidence code
GOs_no_IEA <- GOs[GOs$go_linkage_type != "IEA" & GOs$namespace_1003 == "biological_process",]

dim(GOs_no_IEA)
[1] 111944      5

Then you can continue with goseq.
library(goseq)

All <- read.table(input_All, sep = "\t", header = T)

All_genes <- as.character(All$Geneid)

TopTable <- read.table(input, sep = "\t", header = T)

sig_genes <- as.character(TopTable$Ensembl)

# Make gene vector

gene.vector <- as.integer(All_genes %in% sig_genes)

names(gene.vector) <- All_genes

gene.lengths <- All$Length

pwf <- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
head(pwf)

GO.wall <- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")


text:
Retrieval of genomic position by using biomaRt package
I have a list of several protein names, their primary gene name and their amino acid sequence (as extracted from the 'SEQRES' section in their corresponding PDB files).
I'm looking for the genomic position (triad of translated nucleotides) of a specific amino acid in the sequence (i.e. residue number 64 in a protein sequence of 344 amino acids).
Is it possible to BLAT an amino acid sequence by using biomaRt in R?
Do you suggest other methods to achieve this task?

Do you have the position in the protein from the PDB file? If so, it'd be much easier to use that then the sequence itself.

Yes, the position in the protein and the sequence are extracted from the PDB file. What do you suggest?

Simply use the protein ID and position within the protein, as suggested by Emily_Ensembl. Blasting or blating things would be both slower and more likely to produce incorrect results.

BiomaRt cannot do that I'm afraid. There is an Ensembl REST API endpoint that will get you the genomic location of protein coordinates. It needs an Ensembl peptide ID an input though, so you could use the xref endpoint first to get that. There's a bit of sample code in R on those pages, you could use that to script together if you've got a list of these.


text:
I need to use my code with fasta format. Can anybody help?
I need some help. I have written a program that simulates protein cleavage but I need to modify it so that it produces fasta formatted output.
here is all of my sequence in .fas file
>jgi|Cormi1|6|CCM_00006m.01
MSSFTRAIRPALQGGRRVLLAQRLSTPRIAHPIAATALRSLSTTAPRRSVDLSELQPTPI
THLSEVETAMADAVSKFATEVVLPKAREMDEAEAMDPAVVEQLFEQGLMGIEIPEEYGGA
GMNFTAAIIAIEELARADPSVSVLVDVHNTLCNTALMKWGSDALKKKYLPRLATNTVASF
CLSEPVSGSDAFALATRATETPDGYTITGSKMWITNSMEAGFFIVFANLDPAKGYKGITA
FVVDKDTPGFSIAKKEKKLGIRASSTCVLNFDDVAVPKENLLGKPGEGYKYAIGILNEGR
IGIAAQMTGLALGAFENAVRYVWNDRKQFGTLVGEFQGMQHQIAQSYTEIAAARALVYNA
ARKKEAGEDFVRDAAMAKLYASQVAGRVAGLAVEWMGGMGFVREGLAEKYWRDSKIGAIY
EGTSNIQLNTIAKLLQKEYTA*
>jgi|Cormi1|8|CCM_00008m.01
MLAGSNARILRYIVVAGIILIFFYHFSGSRYDTLPRPVPAPDHAQAASQPSPPNDEATTA
GRPVDSTQKSIPADRDRPNSGFPKSVSASTDAAASTLADIDPGPRANATFVSLARNADVW
DIARSIRQVEDRFNRRFNYDWVFLNDKPFDDNFRNITSALVSGTARFGQIPVEHWSVPDW
IDEDKATKVREDMKERKIIYGDSMSYRHMCRFESGFFFQHPIMKNYDFYWRVEPSIELFC
DIAQDPFRFMKENGKKYSFTISLLEYVDTIPTLWDSAVKFIENYPQHVAADNSMEFLSDD
GGQTYNMCHFWSNFEIGSLEWMRSDAYVDFFNSLDQDGGFFYERWGDAPVHSIAAGLLLN
KTELHFFDEIAYFHMPFTHCPNDVQSRIDLKCHCNPKDNFDWSRYSCTSRYFEMMGMQKP
DGVD*
>jgi|Cormi1|10|CCM_00010m.01
MSTSPSDGTVDVDLPDAPTAPFPLHSDTHTANHVENDEDDDPITASYPVFLNPALPLGRR
LLVLQHPNRTDDHGARAPPSELRLKAESGMVEVDMPLDYNEAYDRDKGLRWGRTLQTSMA
AKNGGSHGLAGGFGFGAVQQRGAGAGAGGARKRNEVDAEDLIDWNEAVRQDKVLRTQTLG
GQYPDADEVQYMVGVFQGNDLHLTPVSSLVHLRPQLHHLDAATQQERQAAAKEAGAAAAS
SSGGATAGPGATGNPAARAIHMSIKTTADGETVTTETIADRLRFVQGEAWRKLRYTDENE
EAAWEVYNESLFLSASEARAPPPGSEKKKKTGSNEEELPYRSLQGAVATFASKWESGQLL
EAVSGMKKPPPVVVVPAAVKKAKAAAEAEAKAKDVEEKRKAAEAQMAATATAAAAAEETR
RPKRAPRGGGGLGTAKRPSRAK*
>jgi|Cormi1|11|CCM_00011m.01
MASATTPILPAEGKRNVLITSALPYSNNLPHLGNVIGSVLSADLYARYCRARGYQTLYVC
GTDQYGTATETKALEEGTDPASLCAKYHALHKSIYDWFRIDFDIFDKTPTQHQTDIVQDI
FTKLWKNGYIYEDETIQPFCPVESHNTFLADRFVEGECSICGDKGARGDQCDACGSILDP
LKPEETGAVDGDAAIEVKATGWLVNPRCKVDGTAPEARKTKHLYLKLDDISEELKAWYKT
TADNGWSTNAKAITQSWIDKGLKPRAISRDLKWGVPIPAVEGLNVEDYARKVFYVWFDAC
IGYISITKTYTDGDNRDGKNWEKWWKNPENVELMQFQGKDNVGFHSIVFPSSLLGTKENW
TKVTRISATEYLNYEQGKFSKSRGIGVFGNNAKDTGIDPDIWRFYLLSRRPELADSEFKW
EEFIDVNNNELLKNVGNFIQRTLKFCVAKMDSTVPAATVSNDVIESHKQSANQILASYIE
HQEANKMRIGISDILTFSALGNKFLQDNGLNNKLLAEEPELCNTVINVAINHVHLLAGLL
YPYMPTTAESIFQQLGQEPNPFIPDVFGFTAVPVGHKIGTPGPLFSVIPAAKLDEWRDAY
GGAELKKQKAIEAEKAAAKKAKKLADKEKKKAKKAETADLKTEGASTAADISKLSLEEKK
*
>jgi|Cormi1|15|CCM_00015m.01
MAEEEVRSTVALVVGALLPALMTVFFNGLPGAVGGSNLPRWGVRLPILLTCAAQLWFMLV
SRVVPEPYLDEIFHIPQAQKYCEGRFFDWDDKITTPPGLYFFSVAVHRIASALRIPQLSV
CDPFSLRAVNFAGVLVVSYLALWCRQALEARQHEAISSPTPARVRAFSQYAVHTAVNIGL
FPLLFFFGGLYYTDVLSTGVVLAAFLNHLRRVGVPRSSVWSDLGTVALGVAALGMRQTNV
FWIVVFMGGLEAVHAVKTLRPARVDQPVMTTLAAQCRYFAWRYSLGDIHDLPVHKAYPDD
MIFTALSIGIAALCNPLRVVRQIWPYLAVLASFAAFVQWNGSVVLGDKSNHVATIHLAQM
LYIWPFFAFFSLPLLLPSGLSFLNMVVSFFRTHARASPTPDAAADGLLRATLRFFHSKVL
WPAYLAATAVLSLAVVRFNTIIHPFTLADNRHYMFYVFRYTIRRGGAVRLALALPYTVCR
WVVWDALGGHNSWSLGNDKPTAFVSHPFWTPPSAKVRRQTSTSYPRPPERASTAESRQAS
EAALKARLEQDSLYLSQQSVSTSTGLIFLLATALSLVTAPLVEPRYFILPWVVWRLMVPA
WRLQDHLRGGILEGSPAAKAVAHWSASYDLRLFLESAWHVVINLATGYVFLYKPYVWRAE
DGTVLEGGNLQRFMW*
>jgi|Cormi1|20|CCM_00020m.01
MTNTIVIVGQVLFPLLTRTTNCIANERRPPSSGGVSGLTSAYLLSKSKGNTVTLVAKHMP
GDYDIEYASPIAGANVLPMNPDADSRWERRTWPEFVWLCRNVPEAGFHFKNALVLRQHKD
QEADPSGGPIDPLFDRDPWYRHLTPDYRELAASELPPGQDSGCSFTSMCLNTPMYLFWLV
GQCLRQGVVFRRAVLHDLTEARRLAYSGRPADVLINATGLGARTLGGVRDAGMLPARGQI
VLVRNDGSSMLNTSAVDDGSGEILYTMQRAAGGGTILGGTYDLGNWESQPCPNVATRIMA
RAVRLDPALAGGKGVAGLDVIRHAVGLRPYRKGGVRVERETLRDEVDGASVHVVHNYGHA
GWGFQGSFGCAEKVVELVNEIRAANGEDLSKEPTLFNFEAAAASKL*
>jgi|Cormi1|21|CCM_00021m.01
MEATNSPARAAREPAVHVRALTFVVDVACGRHALSPLIPVALWLADAVLTCLIIWKIPYT
EIDWVAYMEQITQFVHGERDYTKIKGGTGPLVYPAAHVYTYTGLYYLTGEGENVLLAQQL
FAVLYMVTLGLVMLCYWKAKVPPYIFPLLILSKRLHSVFVLRCFNDCFAVFFLWLAIYSF
QRRLWALGSMAYTWGLGIKMSLLLVLPAVGAVLLHACGFSGALSLAWLMAQVQIVIALPF
IATNTNGYLARAFELSRQFQFEWTVNWRMLGEATFRRREWALLLLGLHAATLALFLATRW
MRPAGKTLRRMVPAWLRGRAALTPDEEIRVSRRVTPDFVLTTMLSANLVGLLFARSLHYQ
FYAYLAWSTPFLLWRALPGGAAVLTYPLWLLQEWAWNVFPSTSESSTAAVAVMATTVLVV
YFGTRNDGLTIKKPAAKKQ*
>jgi|Cormi1|32|CCM_00032m.01
MRLIAAQMLRGTAQRRGFATTPRRLDNYAFVGLGQMGYQMAKNLQSKLQPSDKISLYDIN
TDAMKRLQVEMKAASNGAAVELAASAHDAAKEADTVITVLPEPQHVQGVYKSILTPSLPR
RDRVFIDCSTIDPSTSRAVAASVAAAQQGAFVDAPMSGGVVGATAGTLTFMLGAPASLLP
RLEPALLTMGRKVLHCGDQGAGLSAKLANNYLLALNNIATAEAMNLGIRWGLDPKKLAGV
LNVSTGRCWPSEVNNPVKGVVETAPAGRDYSGGFGIALMKKDLRLAMVAAAEAGADMKLA
DAAFKVYEEASEREDCKGRDFSVVYRYLGGKE*
>jgi|Cormi1|35|CCM_00035m.01
MPGLPTSVDLDECISRLYKKELLAESVIEAICAKTKELLMRESNVVHVQAPVTVVGDIHG
QFYDLIEIFRIGGYCPDTNYLFLGDYVDRGMFSVETISLLVCLKLRYPERVHLIRGNHES
RGVTQSYGFYTECSRKYGNANVWHYFTDMFDFMTLSVVIDNQIFCVHGGLSPSIHSIDQI
KIIDRFREIPHEGPMADLVWSDPDPDRDEFSLSPRGAGYTFGAQVVKKFLAVNNMSHILR
AHQLCQEGFQVLYEDRLSTVWSAPNYCYRCGNMASVLEVNSDGERFFNVFAAAPENDQVK
DMQPGIDKAGDGNALPDYFL*
>jgi|Cormi1|37|CCM_00037m.01
MRYAFAYAALVSTASAHGLVTMIKGDNGVMMPGLGVTDGTPRDCSSNGCGSQADTSIIRN
RDMQQNPLGKTQGTGPVDAAVNVAAFMGTSKNAPKNNGASAGVGQEDDLSGLPAGAAKRH
EERKNMVRGLFGDLLGGAAGGAGAAGGAAGGAAGGLGGLLGGAGGAGGAAGGAAGGLGGL
LGGLGGGGGGTKSNLAAESMVADTAGMGSTQGMPTSDSTGTVSLVYRQINQDGAGPLTAA
VDGTSGGTDASAFQQAEVTKDVPGLGIGGLSLATNTDFPLQVKMPAGMTCEGTVGGVKNV
CIVRVRNSAAAGPFGGSAAFTQNATARKRAIAYRLRKRMEIAREAVPQIDYTKADAAEEL
EDEEDE*
>jgi|Cormi1|65|CCM_00065m.01
MRWHVYWGVVQLVELSNVDKRVAVACHVHTLALPDISLSHHPASSATVTPVIIVRLNRPE
SRHAYTDTMGASLVTAYGLLSADPRVRAIVLTSADASNKFFCAGMDFNVRAQRVEDPLQY
RDSGGTVSIAMHRCQKPIVAAVNGAAVGIGMTMTLPATVRVASTKAKCGFVFAQRGFCME
ACSSFFLPRLIGTAKAMHLVSTGAVYPATHRLFDDLFSEVVEPEDVLPTALKIAEGMAQN
CSIVSTQVMKDMMYSGPKSPEEAHLLESKIFYSLTCPDTDSAEGVQSFLQKRPPKFTGTM
ENQRPFGWPWWSNLDVRPKI*
>jgi|Cormi1|71|CCM_00071m.01
MSRIEKLSISGVRSFSPACREAIQFNTPLTLIVGYNGSGKTTIIECLKYATTGELPPNSK
GGAFIHDPKLCGEKEVLAQVKLQFRSINDRQHVATRSLQLTVKKTTRSQKTLDCSLVVLN
NGEKTTTSTRQAQLDEMIPERLGVSPAILDAVIFCHQDESLWPLSEPSVLKKRFDEIFEA
LKYTKAIENLKVLRKKQVEELGKLQNDEFHNKGNKDRGDRAEKRMTGLQAEIEAAREKCE
AITAHMQETQDKIKQKHEQANSFLQIVQNLGNKREQLDYRQDAVDELRQTIDELPEDDRS
LERQLVQYEESMRRQSEEAQRNKTQYGDLQQDLVKSRKSLSTKLAEQGKHQSDKEKHERQ
LKMRMEMVQEAAQHHGFRGYEDDLTDALIKQFNDKIQKLFAEKKRELERLQKENAAELDK
ASATISEFEGQKAARTQDRVSAKQRQGAIEKRTIVLQNEASSIDVDEGLKAILDGQLGEL
ESKFSQVHREFETAAWDNQISDENAKLWQLESESEKLGRELIECTRLASERAQLDFRRKE
LSDRKLKLTSLTETWRSKLDKVIGGDWSPEDLEAKFQALLAAQNKALQEARQQQEQTRQK
QARVDFKLKTVKETCEQKQTDQIECKKKVLQALQNVRDSAVIEDYTEEVSATEQQIETLR
DDLSLIDALKDYYTKCKRALDAKKTCLLCDRHFDDSQAASIDRLRVKIDRHLNSKTKIEA
EQDLSEASSTLQSLVAVRSQYDTYERLGKELPGLRDERKSIEADHESTERLLEEHDAVVA
AEEDKKRDMNDMSKTVSNITQTLREIKDSEGQVERIASQQQSVGTGRSADEIHDLQGTVS
EQMRASKNRVAKLSNDRQRARDQISTLELEKSELRNKISLTLGQLDKKQDLLGQIQSLKE
EQTHQREAIQRADEELEKVEPRITEARSIRDETLKRGRAKEQVVTESRDSVATSINEIKM
IESDIRDYEDRGGASNLASNQRAIATLEQNITTTETEISELTQRTNKLKQDIDNGDRTKK
NINDNLNYRKHLRMLETLRGEINELEDRNAHEDYERLQGEARALENQSNRLFAERGSVMG
SMKTKDEDLDRLLKEWEMDYKDAKQKYRESHIRVETTKAAIEDLAQCGAAVDKAVMQFHA
MKMAEVNRIAGELWQSTYQGTDIDTILIRSDNEGTTGRRNYNYRLCMVKQDAEMDMRGRC
SAGQKVLASIIVRLALAESFGVNCGLIALDEPTTNLDRDNIKSLAESLHNIIRTRQAQSN
FQLIVITHDEDFLRHMRCSDFCDSFFRVKRDDKQNSVISRESITKIF*
>jgi|Cormi1|75|CCM_00075m.01
MPSSTPLPPPKGFLSYPDVEEIISPRYEPRKDANPVLSGLPLVIASSLLTHSSFLAKIFY
NNAGLQKIKDLPFLDDVPYTFHPLVIPLGDNGPMVELGEAQLNPKNRGDTHHYTASDYHA
MYKSGQVTPLQVTETLLSITSKYAKEASIYQDAWADSHGADQLALDAAKASTERWAAGKA
LGVLDGVPIGVKDDLSVKGYIDHDGMAYKKGVPFFKSAEETVWPVQALQAAGAVVIGKNR
MHELGSDTNGLNPHQGTPTNHLNNQYYPGGSSSGAGSSLGAGVIPIAVGTDAGGSVRFPA
VYNGVYGLKPTHHRTLEMNSTMCVAGPMAANVADLTIAYRIMSQPDPTTSIQSRFAVSQP
PAAGTNKYIGVDRAWWGASDPRVIGACNKALEYFTSQGYELVDITIPYIDEAQTAHGGVT
IAEMTEAARRRGGDKTSWTTLIGPVNKLVLALGLETNAADLLKFNSMRTLIMRHMAWLFQ
KYPGLLILTPTAPFLGWEKKPAHAKYGVSDGDKTFLSMIYVFLANIAGLPALTAPIAYVD
PDQGEGKVCIGMMATGEWGAEEALLAWAADAETYLHEVYPEGRRRPKTWVDVLGLAAKQG
TV*
>jgi|Cormi1|77|CCM_00077m.01
MLSSTITRSRPTATMAKLAGGMLVAGLATKYTIGVAHAESPAAPPRSQVFGRGPALVSLV
LEDAEMVNHNTRRLRFRFQDKEAASGLPLTSALLTYSWPKDSWLPVIRPYTPITTADETG
HLDLLVKQYPNGKSSTHLHSLKPGHTLFFVAPLSGHQWQPNSYKHVTLIAGGAGITPIYQ
LAQGILNNPDDKTAVTVVFGVNSDRDVLFKKEFDEYEQRFPGRFRAVYTVSSPVPGSPHR
KGHVNEALLKEVTPAPGPGNMVFVSGPPPMETALVGSRREPGILEKYGYRKEQVYKF*
>jgi|Cormi1|78|CCM_00078m.01
MAADASFDNLARQFSHLPLIDKLPNCFPDINPVDIYRAHITSILSHITGADKSLIYPVLQ
WTATLDKGDLMLAIPALRIKGKPDDLGKQWVEKWPESPLVEKPVHNGAFISFFFKPEPMA
KTVIPMALQHKQDFGSNPANGLRDQNDPSKGRQKIVIEFSSPNIAKPFHAGHLRSTIIGG
FLANLYRSAGWDVATVNYLGDWGKQYGVLALGFEKYGSEEELVKDPINHLFTVYVKISKE
VAEEKEQADKLKAEGKEEEAQKILDNGVDEKARSYFKKMTEGDEQAIALWKRFRDYSIER
YKQTYARLNIHFDEYSGESKVSEDEMQQAAAKLKEMGLAHEDRGALLMDFTKSVPGKEGK
TLGKAILRKKDGTALYLTRDISELINRHERYKFDKMIYVIASEQDLHVKQFFKIVRMMGR
DDIADKCQHVNFGLVLGMSTRKGTVKFLNDILRDVGDHMHDVMKKNEDKYSQVEDPEQIA
DKLGISSVMVQDMNGKRINNYTFDLSRMTSFEGDTGPYLQYAHARLCSIYRKAELSDEAI
AAADLSLLTEPMAVNLLRVIAQWPDVVQNVLKTLEPTTILTYLFKLTHTLSTSYTHLRIV
GSEPEVMKARLALYNSTRVVLYNGMAILGLTPIERM*
>jgi|Cormi1|79|CCM_00079m.01
MRSRNASPHIGPWRNPARRVRHGTTRNYFDSGCSVWMSSPSRLSAPVPPASPLQKNKFSH
IAIYEQRAATGGVWNTSHVSQEPGFSIPRTTPTSEHEYVVVEDGHDAGARVELVSPVYDN
LETNIPHGLMSYTDLEFPAGTPLFPEHGTVLAYIQQYGRDVEHLVTFETQVRDVRKTTAE
DGRLGWTVETKDLKSGVVSAKSYDAVVAASGHYSDPFVPDIPGIAEFEAVQPGAILHSKF
YRRPEQFTGKKVVVVGNAASGIDISHQIATSAQLPVLISEKDVPGAPPAASNNWSRHVGQ
IVELLPVTRSVRFSSGRIEDGVDAVIFCTGYHYSFPFLGSLAPGITAPNGTYADHLWEHM
LYAPDPTLALLVIPKRIVPFPFAEAQMAVIARIWAGRLNVPNQSERDAWVRRCGAESAPG
ARHTFGYPQDVDYINRLHALSGQARPDAERGLENDGAGKRSPYWGEELSWVRSQIFAIKV
ASRALGERRKECKTLADLGFDFQAAK*
>jgi|Cormi1|87|CCM_00087m.01
MFLLATLFALSALPWTAPAMRSDHLAKLRQETVDMFYHGYNNYMLHAFPEDELRPVSCQP
LTRNPNKPADIGLNDVLGNYSLTLIDSLSTLAILAGGPANEKYTGSQALSDFQDSVAQFV
SYYGDGRQGSSGQGLRSRGFDLDSKVQVFETVIRGVGGLLSAHLFAVGDLPIRGYDPKPD
SKWASEDPLETAPISWHGGFRYDGQLLRLALDLSERLLPAFYTKTGIPYPRVNLRSGIPF
YINSPLHQDSEDTEPVVGNAEITETCSAGAGSLTLEFTVLSRLSGDPRFEQAAKRAFWEV
WERRSDIGLIGNGIDAERGVWIGPHSGIGAGMDSFFEYALKSHILLSGQETPNMTRSHRQ
STTGWLDPNSLHPPLPPEMHSSDAFLEAWHQAHASVKRYIYTDRNHFPYYSNNHRATGQP
YTMWIDSLGAFYPGLLALAGEVDEAVEANLVYTALWTRYGAIPERWSVRESNVEPGIAWW
PGRPEFIESTYHIYRATADPWYLRVGEMVLHDIQRRCWVACGWAGLENVLSGEKQDRMES
FFLGETTKYLYLLFDPEHPLNNLDAAYVFTTEGHPLIIPRRARGSSTPRRRASRTSRQEE
KDVAVYSYYDTRFTNSCPAARTPTEPLAGSNTAARRNLFDVSRFTNLYNTPNIHGPLESF
KTRDEHKGIITQYRATSNYSVFPWTLPPTMLPHNGTCPVPSNRILSWIEFPTVDPAPSLF
ATRGGGGDSNTPLRWYHNVGPTVERVDGLKLQLEQEFNEAAGTDVWMVTHVAGKTVGRHE
NVIFHADHVRHFRDDAFSVLRRKDMVEVVLLVQPEPTPQASRELPPAVMPLPVPGALDEV
SSLDAAPDNTSSTRRPSVRASGQPPATDPDSLFKSILRAVSSALEGSATASTPSPSPLRG
SGSGSGSRDDSSNSDDDVGVASIHSWLANTATGAGAFPLPSVRDALLRNSPDFSAADPAA
SLPWRHVYMAGSACDGPLPEDAPRRHQVIVMRRGGCPFGDKLAHVPNFVPDEASLQLVVV
VDEEGGGHVDDAPRPLLATEQRTPRGTPRLHGIPMVLLPASKGTYELFGEAAAVGLRRKY
VVKSQGYVIENAIVM*
>jgi|Cormi1|88|CCM_00088m.01
MSAGHRVARVHLKEKNTAQDASGRQPSSSSSATAQLASAKVVLGNNLRRFPDFPIPGINF
VDIMPLFANPEAHATLLNALELQITQTFVNTKPDVIVGLDARGFLFGPGLALRLGVPFAA
VRKQGKLPGPCVTAEYVKEYGKDFFQMQEDAIREGQKVLVVDDIIATGGSAKAASDLITQ
LKGHVMGYLFILEIPGLNGRDKLGAAPTTILLQDE*
>jgi|Cormi1|96|CCM_00096m.01
MHLVPKELDKLVITQLGLLAQRRLARGVKLNHAEATALIANNIHELIRDGNHTVADLMAL
GATMLGRRHVIPSVCHTLKEIQVEGTFPTGTYLVTVHNPISTNDGDLARALYGSFLPVPR
DADRLFPAALDPADFAPARQPGAVVAVKGDKITLNADRRRIALRVVSEGDRPIQVGSHYH
FIETNPQLVFDRARAYGFRLDIPAGTSYRFEPGDAKTVTLVEIAGHQVIRGGNGLATGGV
ERWRVDAIVERLQKAGYAHAPEPEVNGVVAGRRRAAAARVPRPYRMDRTAYTVMFGPTTG
DRVRLSSTDLWVKVERDYTTYGDECKFGGGKTLREGMGQASGRPDADSLDLVITNALIID
YTGIVKADIGVKNGIIVGIGKAGNPDVMEGVSPNMIVGSCTDVIAGENKIITAGAIDTHI
HLICPQQAQEALASGVTTFLGGGTGPSAGSNATTCTPGKHLMKQMLQACDALPVNVGITA
KGSDSAPEALREQVAAGACGLKIHEDWGATPSAIDTCLSVCDEMDVQCLIHTDTLNESGF
VESTIASFKNRTIHTYHTEGAGGGHAPDIISVVERANVLPSSTNPTRPFTRNTLDEHLDM
LMVCHHLSKDIPEDVAFAESRIRAETIAAEDVLHDTGAISMMSSDSQAMGRCGEVVLRTW
NTAHKNKLQRGCLPEDDGTGADNFRVKRYVSKYTINPALAQGFSHLVGSVEVGKLADLVV
WDPAWFGTKPHQVIKSGLIAWSQMGDPNASIPTVQPVIGRPMFASFVPAASVVFVSQASV
ASGAVASYGLRKRVEAVRGCRTVGKADMRFNDAMPKMRVDPESYTVEADGEVCTAEPAET
LPLTQSWYGILDNQNAPA*
>jgi|Cormi1|100|CCM_00100m.01
MAARLSAAALAPLRPISSSVCAARQLSQSLLQSQQTRGVAQTYLKKVAEGEERWKARAER
IQNGEQRHVWDILDERGYIKDVAGNTDRIKEVMRIKRIGAYVGIDPTADSLHVGHMLPLM
PLFWLWFHGYPAVTLIGGSTARIGDPTGRLQSREHITNSDISKNITKIHYQLSKLWHNVV
QMRIKYGYEDDWAAKRSLLNNNMWLQGLTLYDFIKRLARDTRIGPMLSRDTAKRKLTDGD
GMSLGEFMYPLLQGWDFWHMYNKIGVQMQIGGSDQYGNIVSGIESLKTIRETEEAPYARM
ESGWQHDPIGFTVPLLTDSAGVKFGKSAGNAIWLDEFKTSAFELYGYFMRRTDDEVEKLL
KLFTFLPMPKIQSIMEEHNADPAKRIAQHNLAFEFVSLIHGSQKALEEAQQHSFRFGGEL
PSIVKEPAASSGIINANTAPRSDIKLPRSVMKLSPAKLLFATGLAPSAAEGQRLVKMQGA
YVAAQPGQKRGLVPGNLNWTPMKMWYPEETSKFLLDDRMLILRKGKHNVRIVELVDDEEW
KASGEFYPGEPYTGEVRRRKQEIIAQAAEDGKTLTMTQVNKMLKKKGRQEGLRVANNPDI
ELPTKTEVRERVRLTHMHKKEERRTRGDEHEW*

My current code:
import re

myprotein = ""
my_protein = []
for protein in myprotein:
    myprotein = protein.rstrip('\n')
    my_protein.append(myprotein)
    my_pro = (''.join(my_protein))

   #cleaves sequence
peptides = re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', my_pro)

print (peptides)

the output is like
MSS
F
TRAIRPAL
QGGRRVLL
AQRL
STPRIAHPIAATA
L
RS
LSTTAPRRSVD
L
SE
L
QPTPITH
L
SEVETAMADAVSK
F
ATEVVL
PKAREMDEAEAMDPAVVEQ
L
F
EQG
L
MGIEIPEE
Y
GGAGMN
F
TAAIIAIEE
L
ARADPSVSV
L
VDVHNTL
CNTA
L
MK
W
GSDA
L
KKKYLPRL
ATNTVAS
F
C
L
SEPVSGSDA
F
A
L
ATRATETPDG
Y
TITGSKM
WITNSMEAG
F
F
IV
F
AN
L
DPAKG
YKGITA
F
VVDKDTPGF
SIAKKEKK
LGIRASSTCV
L
N
F
DDVAVPKENL
L
GKPGEG
Y
K
Y
AIGI
L
NEGRIGIAAQMTG
L
A
L
GA
F
ENAVRY
VW
NDRKQFGT
L
VGE
F
QGMQHQIAQS
Y
TEIAAARA
LV
Y
NAARKKEAGED
F
VRDAAMAK
L
YASQVAGRVAG
L
AVE
W
MGGMG
F
VREGL
AEK
Y
WRDSKIGAI
Y
EGTSNIQ
L
NTIAK
L
LQKE
YTA

I need every sequence on a separate line, starting with the sequence header like:
>jgi|Cormi1|6|CCM_00006m.01 frag.1
MSS
>jgi|Cormi1|6|CCM_00006m.01 frag.2
F
>jgi|Cormi1|6|CCM_00006m.01 frag.3
TRAIRPAL

How could I make this work?
import sys
import re

sequence = open(input('Enter input filename: '),'r')
def digest(sequence):
    return re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', sequence)

seq = ""
header = None
for line in sys.stdin:
    if not line.startswith('>'):
        seq = seq + line.rstrip()
    else:
        if header:
            for idx, fragment in enumerate(digest(seq).split('\n')):
                sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))
            seq = ""
        header = line.rstrip()

# final pass
for idx, fragment in enumerate(digest(seq).split('\n')):
    sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))

Use BioPython for this with SeqIO.write

what do you mean? using for which step?

You can use it for both reading in amd writing out your sequences in a preformatted FASTA format. There is rarely a good reason to write your own fasta parser.

One way is to store the header as a separate string from your protein residues. Use a list comprehension to build a list of formatted header-fragment pairs, and then print them to standard output.
#!/usr/bin/env python

import sys
import re

myprotein = ">jgi|Cormi1|6|CCM_00006m.01\nMSSFTRAIRPALQGGRRVLLAQRLSTPRIAHPIAATALRSLSTTAPRRSVDLSELQPTPI\nTHLSEVETAMADAVSKFATEVVLPKAREMDEAEAMDPAVVEQLFEQGLMGIEIPEEYGGA\nGMNFTAAIIAIEELARADPSVSVLVDVHNTLCNTALMKWGSDALKKKYLPRLATNTVASF\nCLSEPVSGSDAFALATRATETPDGYTITGSKMWITNSMEAGFFIVFANLDPAKGYKGITA\nFVVDKDTPGFSIAKKEKKLGIRASSTCVLNFDDVAVPKENLLGKPGEGYKYAIGILNEGR\nIGIAAQMTGLALGAFENAVRYVWNDRKQFGTLVGEFQGMQHQIAQSYTEIAAARALVYNA\n"

# parse record
my_protein = []
for line in myprotein.split('\n'):
    if not line.startswith('>'):
        my_protein.append(line)
    else:
        header = line

# build sequence
my_pro = ''.join(my_protein)

# cleave sequence
peptides = re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', my_pro)

# build a formatted list
records = ["%s frag.%d\n%s\n" % (header, idx + 1, fragment) for idx, fragment in enumerate(peptides.split('\n'))]

# print list
for record in records:
    sys.stdout.write("%s" % (record))

This gives:
>jgi|Cormi1|6|CCM_00006m.01 frag.1
MSS
>jgi|Cormi1|6|CCM_00006m.01 frag.2
F
>jgi|Cormi1|6|CCM_00006m.01 frag.3
TRAIRPAL
...
>jgi|Cormi1|6|CCM_00006m.01 frag.89
LV
>jgi|Cormi1|6|CCM_00006m.01 frag.90
Y
>jgi|Cormi1|6|CCM_00006m.01 frag.91
NA

If you have multiple records in one file called seqs.fa, then you could use standard input to this script (e.g., named digest.py) to process them all at once:
#!/usr/bin/env python

import sys
import re

def digest(sequence):
    return re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', sequence)

seq = ""
header = None
for line in sys.stdin:
    if not line.startswith('>'):
        seq = seq + line.rstrip()
    else:
        if header:
            for idx, fragment in enumerate(digest(seq).split('\n')):
                sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))
            seq = ""
        header = line.rstrip()

# final pass
for idx, fragment in enumerate(digest(seq).split('\n')):
    sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))

Creating a function called digest lets you change the digestion rules in one place. It reduces code bloat and also reduces errors from changing things in one place and not another.
You could use this script like so:
$ ./digest.py < seqs.fa > seqs.digested.fa

The file seqs.digested.fa would look like:
$ cat seqs.digested.fa
>jgi|Cormi1|6|CCM_00006m.01 frag.1
MSS
>jgi|Cormi1|6|CCM_00006m.01 frag.2
F
>jgi|Cormi1|6|CCM_00006m.01 frag.3
TRAIRPAL
...
>jgi|Cormi1|100|CCM_00100m.01 frag.143
PTKTEVRERVRL
>jgi|Cormi1|100|CCM_00100m.01 frag.144
THMHKKEERRTRGDEHE
>jgi|Cormi1|100|CCM_00100m.01 frag.145
W*

Thanks Alex Reynolds but  I need more function to my program that make it can read file (.fas) and run overall sequence in file (each sequence has different name on header)

See edited answer.

sorry for bother you again but I don't get how to make this work. So after creating the function "digest" where should I put the script to? It would be great if you explain step by step.( I am new for this)

Copy the script to a file called digest.py. Make the script executable via chmod +x digest.py. Say you have a file containing your sequences at seqs.fa. Run ./digest.py < seqs.fa > seqs.digested.fa to make a file in the format you want.

so I just upload the picture of the way I think you told me. Is it right or not sorry for my dumbness I really don't know how to make this work. :'(

I haven't used Windows in a long time, but I wish I could help! Is there a Windows-based Python tutorial that might help? You might also try running things through Cygwin.

https://docs.python.org/3.4/using/windows.html

ummm is there anyway to finish this only use in python? I'm thinking about open() or something like that.

If you use Windows PowerShell, maybe you could do something like Get-Content seqs.fa | ./digest.py > seqs.digested.fa. See: https://stackoverflow.com/questions/11447598/redirecting-standard-input-output-in-windows-powershell

seem like I can't done with cmd or powershell or I don't have skill for that haha so I think the best way is using something just in python so I put the line: sequence = open(input('Enter input filename: '),'r') above def digest like in edited. so how can I make this work by this way?

If you’re on windows 10, install the Bash on Ubuntu application that was released as part of the anniversary update a few months ago. It essentially gives you a full Unix virtual machine on windows, and integrates quite nicely. Save yourself a lot of headaches in the future and switch to unix based systems ASAP.


text:
Which data format for stranded coverage data
I need tools I am writing to output data in terms of score for each base in the genome. The format needs to be indexed for fast random access.
The obvious choices are bigWig and tabix-index bedgraph (indexed bedGraph seems quicker for computation, bigWig better for visualization), but the problem is that I need the data to be stranded - i.e. I need to record which strand the data belongs to (and perhaps even have a score on the positive and negative strands).
So my choices are:

Invent a new format
Bodge one of the existing formats (this is almost as bad as inventing a new format, if not worse). For example: bedGraph, except with two score columns - one for + and one for -, or have a bedGraph with with a strand column. This doesn't seem ideal because other downstream tools won't be able to process it.
Have seperate files for + and - strand scores. This is unweidly and not great if you sometimes don't care about strand, or want to visualise on a single track.
Use a magic format I've not heard of but that someone is going to point out in the answers...

None of these are ideal. As they either make tool output files incompatible with other tools, or don't encode all the information I want to. Is there a an established format I can use for this?

Since you lose many benefits of creating a new format: I found that having three files — forward, reverse, unstranded/merged — works best.

So far I've been working with three files: a full bed6 with a line per non-zero coverage base. This I index with tabix and use internally. Two bigWigs - forward and reverse, which I use for external programs. The difference in performance is quite impressive. On one benchmark what makes about 18000 queries, querying BAM and processing to stranded depth takes 45 seconds, querying bigWig takes 90 seconds and querying the indexed bed takes 25 seconds.


text:
How to convert the given mathematical computation (on biological problem) to mathematical fomula, equation?
I have crossposted this question in maths StackExchange. The problem is dominantly mathematical (this question) but the application of the problem is mainly biological. Hoping that people in this forum have faced similar problems, I am posting this in bioinformatics forum to get some ideas and probably solutions to it.
I am working on a problem that allows haplotype phasing. I have developed this computation method (expressed in detail below) and have developed a python code (not shown here) to solve the issue.
Now, I am having a hard time trying to translate this computation in mathematical and/or statistical language. I have revisited maths and I am able to comprehend the process at several places, but I cannot translate those ideas onto my computation (or algorithm).

Here is my computation (step by step) with a very workable example. I have tried to make it as comprehensive as possible, let me know if there is any confustion:
Fig 01: Example of the haplotype file
contig   pos     ref   all-alleles   ms01e_PI   ms01e_PG_al   ms02g_PI   ms02g_PG_al   ms03g_PI   ms03g_PG_al   ms04h_PI   ms04h_PG_al   ms05h_PI   ms05h_PG_al   ms06h_PI   ms06h_PG_al
2      15881764   .      .           4           C|T           6           C|T           7           T|T           7           T|T           7           C|T           7           C|T
2      15881767   .      .           4           C|C           6           T|C           7           C|C           7           C|C           7           T|C           7           C|C
2      15881989   .      .           4           C|C           6           A|C           7           C|C           7           C|C           7           A|T           7           A|C
2      15882091   .      .           4           G|T           6           G|T           7           T|A           7           A|A           7           A|T           7           A|C
2      15882451   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|A
2      15882454   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|T
2      15882493   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|T
2      15882505   .      .           4           A|T           4           T|A           7           T|T           7           T|T           7           A|C           7           A|T

This (fig 01) is my main data file. The main idea is to take a pool of several samples which have short phased haplotype blocks represented as PI (phased index) and  PG_al  (phased genotype). For any sample, the site where two consecutive haplotype blocks are not joined represents the break point (and has different PI values). In the above example only sample ms02g has break point and needs to be phased. All other samples have haplotype block that bridges these two consecutive haplotype blocks and contains the data required to extend the phase state of sample ms02g.
Fig 02: Representing a break-point in the sample: ms02g
contig    pos    ms02g_PI    ms02g_PG_al
2     15881764     6         C|T
2     15881767     6         T|C
2     15881989     6         A|C
2     15882091     6         G|T
                   ×——————————×—————> Break Point
2     15882451     4         T|C
2     15882454     4         T|C
2     15882493     4         T|C
2     15882505     4         T|A

So, in the above haplotype file there is a breakpoint in sample ms02g at the position 15882091-15882451. For sample  ms02g(PI-6) the haplotypes are  C-T-A-G and  T-C-C-T. Similarly, at PI-4 the haplotypes are T-T-T-T and C-C-C-A. Since the haplotype is broken in two levels, we don’t know which phase from level-6 goes with which phase of level-4. But, all other samples have full haplotype intact that bridges this position. We can therefore use this information from other samples to join the two consecutive haplotype in sample ms02g.
Using markov-chain transition to extend the haplotype block:
Since, all other samples are completely phased bridging that breakpoint, I can run a markov-chain transition probabilities to solve the phase state in the sample ms02g . To the human eye/mind you can clearly see and say that left part of ms02g (PI-6, i.e C-T-A-G is more likely to go with right block of PI-4 C-C-C-A), thereby creating the extended haplotype block as  C-T-A-G-C-C-C-A and T-C-C-T-T-T-T-T.
Below, I show how I can apply the first order markov transition matrix  to compute the likelyhood estimates, calculate the log2Odds and then assign and extend the haplotype in proper configuration. And, I feed this logic to the computer using python.
Calculation of likelyhood estimates using markov transition:
Step 01:   prepare required haplotype configuration

The top PI-6 is Block 01 and the bottom PI-4 is Block 02.
The phased haplotype in the left within each block is Hap-A and on the right is Hap-B.

Fig 03: representation of the haplotype breakpoint and block assignment
ms02g_PI    ms02g_PG_al
  6         C|T \
  6         T|C | Block 1
  6         A|C |
  6         G|T /
  ×——————————×—————> Break Point
  4         T|C \
  4         T|C |
  4         T|C | Block 2
  4         T|A /
           ↓   ↓
       Hap-A   Hap-B

So, the two consecutive blocks can be extended in one of the two possible haplotype configurations:
Parallel Configuration:
 Block01-HapA with Block02-HapA, so B01-HapB with B02-HapB
Vs. Alternate Configuration:
 Block01-HapA with Block02-HapB, so B01-HapB with B02-HapA

Step 02: Compute transition matrix and estimate likelihood for each configuration.
Fig 04: Representation of the allele transition matrix (from alleles of former block-01 to alleles of later block-02).
Possible
transitions     ms02g_PG_al
 │              ┌┬┬┬  C|T \
 │              ││││  T|C | Block 1
 │              ││││  A|C |
 │              ││││  G|T /
 └────────────> ││││   ×—————> Break Point
                │││└> T|C \
                ││└─> T|C | Block 2
                │└──> T|C |
                └───> T|A /
                     ↓   ↓
                 Hap-A   Hap-B

**note: In this example I am showing transition only from first
  nucleotide of Block01 to first nucleotide of Block02. Actually we
  prepare transition from all nucleotides of B01 to all nucleotide of
  B02 in both parallel and alternate configuration.

I count the number of transitions from each nucleotide of PI-6 to each nucleotide of PI-4 for each haplotype configuration across all the samples and convert them to transition probabilities.
And multiply the transition probabilities from the first nucleotide in PI-6 to all nucleotides of PI-4. Then similarly multiply the transition probability from 2nd nucleotide of PI-6 to all nucleotides in PI-4, and so on.
When transition probabilities are calculated for all possible combination (from each position of PI-6 to each position of PI-4), then I compute the cumulative transition probabilities for each possible haplotype configuration.

Fig 05 : Representation of nucleotide counts (emission counts) at positions 15881764 and 15882451.
pos\allele     A    T    G    C
15881764       0    8    0    4
15882451       1    7    0    4

Fig 06 : Representation of transition matrix counts (from pos 15881764) to (pos 15882451).
This transition matrix is computed from nucleotides (A,T,G,C) at block01 to nucleotides (A,T,G,C) at block02 for all the positions. Transition counts are then converted to transition probabilities.
from     to
          A    T    G    C
  A       0    0    0    0
  T       1    6.5  0    0.5
  G       0    0    0    0
  C       0    0.5  0    3.5

Note: if the PI matches between two blocks the transition are counted
      as 1, else 0.5. Sample ms02g itself is also taken as an observation.
Step 03: Compute the maximul likelihood for each configuration
Fig 06 : Likelihood estimate for parallel configuration using transition counts (probabilities)
Parallel configuration:
Block-1-Hap-A (C-T-A-G) with Block-2-Hap-A (T-T-T-T)
  CtT × CtT × CtT × CtT = (0.5/4)*(0.5/4)*(0.5/4)*(0.5/4) = 0.000244
+ TtT × TtT × TtT × TtT = (0.5/2)*(0.5/2)*(0.5/2)*(0.5/2) = 0.003906
+ AtT × AtT × AtT × AtT = (0.5/3)*(0.5/3)*(0.5/3)*(0.5/3) = 0.0007716
+ GtT × GtT × GtT × GtT = (0.5/2)*(0.5/2)*(0.5/2)*(0.5/2) = 0.003906
——————— ————————— ——————— ————————— Max Sum (likelihoods) = 0.008828
                                    Average (likelihoods) = 0.002207

Block-1-Hap-B (T-C-C-T) with Block-2-Hap-B (C-C-C-A)
  TtC × TtC × TtC × TtA = (0.5/8)*(0.5/8)*(0.5/8)*(0.5/8) = 0.00001526
+ CtC × CtC × CtC × CtA = (2.5/10)*(2.5/10)*(2.5/10)*(2.5/10) = 0.003906
+ CtC × CtC × CtC × CtA = (1.5/8)*(1.5/8)*(1.5/8)*(1.5/8) = 0.001236
+ TtC × TtC × TtC × TtA = (0.5/4)*(0.5/4)*(0.5/4)*(0.5/4) = 0.000244
——————— ————————— ——————— ————————— Max Sum (likelihoods) = 0.0054016
                                    Average (likelihoods) = 0.0013504

note:
      - "AtC" -> represent "A" to "C" transition
      - "+" represents the summation of the likelyhoods
Fig 07 : Likelihood estimate for alternate configuration using transition counts (probabilities)
Alternate configuration:
Block-1-Hap-A (C-T-A-G) with Block-2-Hap-B (C-C-C-A)
  CtC × CtC × CtC × CtA = (3.5/4)*(3.5/4)*(3.5/4)*(3.5/4) = 0.5861
+ TtC × TtC × TtC × TtA = (1.5/2)*(1.5/2)*(1.5/2)*(1.5/2) = 0.3164
+ AtC × AtC × AtC × AtA = (2.5/3)*(2.5/3)*(2.5/3)*(2.5/3) = 0.4823
+ GtC × GtC × GtC × GtA = (1.5/2)*(1.5/2)*(1.5/2)*(1.5/2) = 0.3164
——————— ————————— ——————— ————————— Max Sum (likelyhoods) = 1.7012
                                    Average (likelihoods) = 0.425311

Block-1-Hap-B (T-C-C-T) with Block-2-Hap-A (T-T-T-T)
  TtC × TtC × TtC × TtA = (6.5/8)*(7.5/8)*(7.5/8)*(6.5/8) = 0.5802
+ CtC × CtC × CtC × CtA = (6.5/10)*(7.5/10)*(7.5/10)*(6.5/10) = 0.237
+ CtC × CtC × CtC × CtA = (5.5/8)*(6.5/8)*(6.5/8)*(6.5/8) = 0.36875
+ TtC × TtC × TtC × TtA = (3.5/4)*(3.5/4)*(3.5/4)*(2.5/4) = 0.4187
——————— ————————— ——————— ————————— Max Sum (likelyhoods) = 1.60465
                                    Average (likelihoods) = 0.4011625

note:
  - the sum of the likelihoods > 1, in the above example.
  - So, we can rather use the product of the likelyhoods.

Fig 08 : Likelhood estimate of Parallel vs. Alternate configuration.
Likelihood of Parallel vs. Alternate configuration

= likelihood of Parallel config / likelihood of Alternate config

= (0.002207 + 0.0013504)/ (0.425311 + 0.4011625)
= 0.0043043  (i.e 1/232)
Therefore, haplotype is 232 times more likely to be phased
in "alternate configuration".

log2Odds = log2(0.0043043) = -7.860006
So, with our default cutoff threshold of "|5|", we will be extending
the haplotype in "alternate configuration"

Final Output data:
contig   pos   ref   all-alleles   ms02g_PI   ms02g_PG_al
2   15881764   .    .                  6      C|T
2   15881767   .    .                  6      T|C
2   15881989   .    .                  6      A|C
2   15882091   .    .                  6      G|T
2   15882451   .    .                  6      C|T
2   15882454   .    .                  6      C|T
2   15882493   .    .                  6      C|T
2   15882505   .    .                  6      A|T

So, how do I translate my overall computation process in mathematical and statistical language.
I am looking to represent computation of:

emission counts (and probabilities)
transition counts (and probabilities)
summation of the transition counts for each configuration
computation of likelihood ratio and it's conversion to log2Odds

I am trying to represent my computation (algorithm) in the way it's shown in these links, but have a hard time doing it. I've read, this, and this, and this.
FYI - I am biologist but have managed to learn python and write a python program to work it out. The code isn't relevant in this question though.
Now, it's my time to work out some maths and represent my algorithm.

Hi, I have edited a bit the question, I think that now is easier to read, feel free to edit further or roll back if you disagree. what is encoded in the haplotype columns (Fig 01) ? It seems that you have like several samples in the same file, and you are only interested in one sample. Have you looked at how do the mathematics represent markov chain  models?

@Llopis: Yes, I am trying to solve one sample at a time; in this case ms02g - because this is the one that has a breakpoint (i.e different PI blocks). Other samples have full length haplotype, so they actually serve has a observation data for preparing the markov transition matrix.

@Llopis: Yes, I looked into how to represent markov model. I can comprehend at several places but all. So, I cannot translate. Also, my data has PI values too which needs to be considered while computing transition counts (and probabilities).

Are you asking how to take what you want to do and write a program to do it, or are you asking, given an algorithm, how to do a write-up of your algorithm to put in a paper?

I have already written a program. I want to be able to write it in a paper. I have left maths for more than a decade and seems like I will be quite slower if I tried to re-read the whole thing again. So, any way you can help is appreciable.

I am not sure I understood all, but what you do is like doing a Markov chain but instead of saying that the previous position(s) decides the fate of the next one(s) is the block what determines which other block follows it.
Mathematically I think that could be explained in terms of linkage disequilibrium (and wikipedia page), but instead of using a single position you use a block. I am not familiarly enough with the field to see if a similar (or this method) has been previously described, but I think you could find inspiration about how to describe it.
A related field of work that also deals with this kind of data (and problems) is structural chromatin, where they work with TADs (their blocks of Topologically Associating Domains) and study which ones are near (spatially instead of probabilistic) with others regions, that field could be also interesting to look for formulas and how do they present this kind of analysis.

Thanks. But, my problem is entirely different. I have looked into markov chain. That is how I was able to work out that algorithm and write a program. Now, I need a mathematical way to represent my algorithm. I think I can go ahead and read TAD paper in detail - for which I don't really have much time. I am sure there is a better way on working this out. Thanks

Then could you explain how your problem is different to linkage disequilibrium? When you say you "need a mathematical way to represent my algorithm" do you refer to a mathematical formula like $y = a +b*x$?

My problem is a more complicated than LD, though the very first principle of LD does apply in some sense. In LD you only have 2 genes, 2 alleles and trying to see the segregation pattern. And, even Markov Chains aren't quite involved. It is very easy to represent LD mathematically. On contrary my method has two blocks (length "n" and "m") and then two haplotypes within each block. Then the transitions counts from each "n" to each "m" which are summed. Yes. I just want a mathematical representation. Look for the links in my main question I've read, this, and this, and this......


text:
counting individuals with the present gene and absent gene for each population
I have a folder containing 500 files from different individuals:
$ ls
HI.3577.005.Index_25.QPG4_KS2W-35_all_PA
HI.3577.005.Index_27.QPH4_SD1W-13_all_PA
HI.3577.005.Index_9.QPA2_KS2W-8_all_PA
HI.3577.007.Index_1.Q004B_SD1A-1_all_PA
HI.3577.007.Index_21.Q076C_SD2W-23_all_PA
HI.3577.007.Index_23.Q079_KS1A-15_all_PA

Each of these individuals belongs to a population (in the example above individuals belong to KS2W, SD1W, KS2W, SD1A, SD2W, and KS1A respectively). Each of my files looks as follows (this sample belongs to from "MK4" population):
$ head HI.3579.001.Index_9.Q089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 0
Ha16_00023451 0
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

or another sample from MK4 population
$ head HI.221.009.Index_9.P089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 1
Ha16_00023451 1
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

the first column is genes and the second columns refers to the presence of that gene (coded as 1) or absence of gene (coded as 0).
For all individuals belonging to each of the populations, I would like to count the total number individuals that have the gene and individuals that do not have the gene.
For example, for the 2 examples above, as they both belong to MK4 population my desired output is going to looke like this:
$ head MK4
Ha8_00040788   2  0
Ha4_00024045   2  0
Ha4_00025366   2  0
Ha16_00022130  1  1
Ha16_00023451  1  1
Ha8_00040789   2  0
Ha4_00025367   2  0
Ha4_00024046   0  2

The first column is the gene, the second column is the number of individuals in the MK4 population that have the gene and the third column is the number of individuals that do not have that gene.
I have made a file containing population IDs
$ head pop.txt
KS2W
SD1W
SD1A
KS1A

Can anyone suggest me how I could do this? thanks

It's unlikely that genes are absent in individuals; do you mean SNPs?

Why not use R??

@b.nota Some people have programming-language preferences.

I can show you how I would approach it in R. You'll have to adjust and customize it to your own data of course. I only had two files to play with, from your example. Especially for the merge step you should build a for loop around it so that it will merge every next file automatically, since you have 500 files.
Here an example for the KS2W files:
myFiles <- dir(pattern="KS2W")

myFiles
[1] "HI.3577.005.Index_25.QPG4_KS2W-35_all_PA.txt"
[2] "HI.3577.005.Index_9.QPA2_KS2W-8_all_PA.txt"

fileNames <- list()

for(i in 1:2){
    fileNames[[i]] <- read.table(myFiles[i], sep = " ")
}

merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

present <- rowSums(merged[,2:3])

zeros <- rowSums(merged[,2:3] == 0)

output <- cbind(present, zeros)

rownames(output) <- merged[,1]

output
              present zeros
Ha8_00040788        2     0
Ha4_00024045        2     0
Ha4_00025366        2     0
Ha16_00022130       1     1
Ha16_00023451       1     1
Ha8_00040789        2     0
Ha4_00025367        2     0
Ha4_00024046        0     2
Ha16_00022131       2     0
Ha16_00023452       2     0

To add more files with merge, you could add something like this:
merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

for(i in 3:30){
  merged <- merge(merged, fileNames[[i]], by.x = "V1", by.y = "V1", all = F, sort = F)
}

You'll have to adjust all the numbers and coordinates to fit your data.

The following python script takes a population name as argument, searches for all files for that population (assuming some regularity in the file names), and records the counts of "0" and "1" for all the genes encountered:
#!/usr/bin/env python3

from sys import argv
from glob import glob
from collections import defaultdict

filenames = glob("HI*Index_*_{pop}-*_all_PA".format(pop=argv[1]))

gene_counts = defaultdict(lambda : defaultdict(int))

for filename in filenames:
    with open(filename, "r") as genes_file:
        for line in genes_file:
            (gene, num) = line.strip().split()
            gene_counts[gene][num] += 1

for (gene, counts) in gene_counts.items():
    print(gene, counts["1"], counts["0"], sep="\t")

With the two example files you give being present in the current directory:
$ ./count_gene.py MK4
Ha8_00040788    2   0
Ha4_00024045    2   0
Ha4_00025366    2   0
Ha16_00022130   1   1
Ha16_00023451   1   1
Ha8_00040789    2   0
Ha4_00025367    2   0
Ha4_00024046    0   2
Ha16_00022131   2   0
Ha16_00023452   2   0

I haven't tested with python 2. If you use python 2, you should at least add from __future__ import print_function as first import, or adapt the print part.
Then, to run on all populations:
for pop in $(cat pop.txt)
do
    ./count_gene.py ${pop} > ${pop}
done


text:
Find specific functional genes from a list?
I have compiled a list of genes that are distinct from 2 different conditions.
From this list (of about 600 human genes), I want to find out if any of the genes in this list have acetyl transferase function i.e. are acetylases.
Is there an online platform where I could say put in this list of genes and the output I can select genes which have this property/function?
I.e some sort of gene annotation for functions for a list of human genes?

Have you looked into any services yet?  Did they not met your needs?

BioMart. Filter by your list of genes and a GO term linked to acetyl transferase activity.


text:
Cufflinks Error: sort order of reads in BAMs must be the same
I am running Cufflinks for transcriptome assembly using the .bam file generated by Hisat2. I tried both bam and sorted bam files
cufflinks --no-update-check -o ./ -p 15 accepted_hits.bam

it gives the following error:
Error: sort order of reads in BAMs must be the same

UPDATE: It is working fine now. It seems there were some issues with the reference file. The BAM files generated from Tophat2 also gave the same error. I downloaded the updated reference file and it worked great.

Maybe post some code so that we can follow what you're doing?

Hello again, could you please add also the Hisat2 command? We are really trying to help, but to do that, we need more info. (do not hesitate to write there everything you tried)

Instead of an update the preferred method would be to create an answer and then accept it.

Just add --dta-cufflinks to your Hisat2 command, so that the output alignment file provides the attributes needed by Cufflinks (XS flags). This should do the job.
From manual:

Report alignments tailored specifically for Cufflinks. In addition to
  what HISAT2 does with the above option (--dta), With this option,
  HISAT2 looks for novel splice sites with three signals (GT/AG, GC/AG,
  AT/AC), but all user-provided splice sites are used irrespective of
  their signals. HISAT2 produces an optional field, XS:A:[+-], for every
  spliced alignment.

If it was what you were looking for, please mark the answer as 'accepted' too. Future googlers will thank you ;)

Sure, I will do that. I am still running hisat2 with that option.

It is still not working !!

I need more detail. It worked for me. Please then edit the question wit exact command line, what you tried and possibly a minimal example to be run. Also, be sure you deleted previously generated files/directories from Cufflinks and Hisat2 before relaunching the clean pipeline.

Btw by 'more details' I meant also to paste in the question the command you used to generate the alignment file with Hisat2, also including the versions of both Hisat2 and Cufflinks that you are using.

Accepted your answer :) I just updated the original question with how the issue got resolved.


text:
How does htslib/samtools access optional BAM fields?
I am confused regarding how various software libraries deal with optional fields in a BAM:
Based upon the BAM specification, there are 11 mandatory fields to a BAM:
QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL

The BAM files of today normally have many optional fields as well, which must be of the format TAG:TYPE:VALUE, whereby TAG is a two-character string that matches /[A-Za-z][A-Za-z0-9]/ and TYPE is a single case-sensitive letter which defines the format of VALUE. Each TAG can only appear once in one alignment line.
Certain BAMs may have no optional fields, while others have several.
My understanding is that optional tags may be in a different order from row to row, and then sometimes they simply do not exist.
My questions:

How does htslib/samtools (e.g. samtools view) access these optional fields, which differ in number from BAM to BAM (possible row to row)?
In htslib/sam.h, a BAM alignment struct is defined as follows:
/*! @typedef
 @abstract Structure for one alignment.
 @field  core       core information about the alignment
 @field  l_data     current length of bam1_t::data
 @field  m_data     maximum length of bam1_t::data
 @field  data       all variable-length data, concatenated; structure: qname-cigar-seq-qual-aux
typedef struct {
    bam1_core_t core;
    int l_data;
    uint32_t m_data;
    uint8_t *data;
#ifndef BAM_NO_ID
    uint64_t id;
#endif
} bam1_t;

and followed by several methods which access mandatory BAM fields, e.g.
/*! @function
 @abstract  Get the name of the query
 @param  b  pointer to an alignment
 @return    pointer to the name string, null terminated
 */
#define bam_get_qname(b) ((char*)(b)->data)

How does htslib access each of the optional fields, given this structure? Does it keep the order of the optional fields, or is this information unavailable?
When a user runs htslib or even samtools view, how does the BAM row "end", i.e. what data structure exists to show that the BAM row has no more optional fields, and that the algorithm should "move on" to the next BAM row.

Any clarification for how the source code works in this case is much appreciated.

Cross-posted on Reddit

How does htslib access each of the optional fields, given this structure?

htslib uses a macro, bam_get_aux, for this purpose. You can check the source code of this macro — it essentially just advances the data pointer of the current record to the start of the optional data, which is after all other data. But you shouldn’t care about this, and just use the relevant exported functions for reading optional (“auxiliary”) data:

bam_aux_get returns a pointer to an optional record, given by its two-character tag name.
bam_aux2‹C› converts the pointer returned by the previous function to the actual type of the optional tag; here, ‹C› is a type identifier: i for an integer, Z for a zero-terminated character string, etc. Refer to the documentation and the htslib/sam.h file for a complete list.

Does it keep the order of the optional fields, or is this information unavailable?

As long as you don’t modify optional fields, the order is kept, but unspecified. Do not rely on the order!

When a user runs htslib or even samtools view, how does the BAM row "end"

BAM records store a header that contains the record size.

@EB2127 No, it won’t necessarily be uniform. And if there are algorithms that depend on the order, then these algorithms are broken because they are breaking the BAM specification. As such, there is no way to set the order. If you need ordered tags (but what for?!) then you need to store the order inside the tag data itself (e.g. by prefixing the tag data with ordinals).

Thanks, I appreciate the help. Another question: "You can check the source code of this macro — it essentially just advances the data pointer of the current record to the start of the optional data, which is after all other data.". Could you flesh out the details for this? This may help me better understand how these macros work. e.g. bam_get_qname()  looks like it's accessing the 'data' member of  struct bam1_t, which is of the form "qname-cigar-seq-qual-aux". Then 'bam_get_cigar()' is advancing the data pointer of the current record to the length of the query name, core.l_qname?

In a BAM file, all of the auxiliary tags (the optional fields you refer to) are held what is effectively a void * field, one after the other. Each entry within this field is structured as a two chars storing the tag name (e.g., NM, followed by a char storing the type (e.g., Z or i), followed by the data in whatever format is appropriate (char, float, etc). Strings in the aux tags are NULL terminated, so their lengths aren't explicitly stored.
To access a given field, htslib has to iterate over all of them, in whatever order they happen to be stored, until it finds the correct one. The output order with samtools view is the order in which they're stored.
The size of *data is stored as l_data, so once it hits that then it knows that the record is finished.

Regarding point 1, this is a standard C trick of casting between data types. So, while *data is actually a uint8_t * (aka, a byte array), the first two bytes of each tag are chars and whatever is after them (well, after the byte specifying the type) is then cast to the appropriate type (be it float or char* or what not). This is all probably clearer if you look at the code.

Thanks for the help; this was very useful. Could you detail #1 a bit more? That is, how the optional fields are stored with a generic pointer type?

@EB2127 I've added a bit more detail as requested. It's best to just look at the source code at this point though. Ignore the stuff on the B type, since that will likely just be confusing at first.


text:
Error while Importing fastq sequences to jMHC software
I am trying to analyze my .fastq files (MiSeq exon reads) with the jMHC software, but output files I have obtained have no results.
I think I might be doing something wrong during the first step: "connect to database".
Because on the "import/extracting" step, I open all my .fastq files, click "Start", and it shows this message:

"Success! Processed: 139, Imported: 0".

So I assume it is reading the files OK, but not importing them to the "database"? Could it be the extension of my files? The software says it operates on FASTA files. Are .fastq files not compatible?
Anybody that has successfully run this program could give me a hand?

Yes, use .fasta , you can try to convert your fastq files with NextGen Workbench for example.

Hello, thanks for the answer. It would be useful if you would provide link to the tool you suggest.


text:
Which R package to use for differential analysis with TPM values?
I'm using hisat2, stringtie tools for the RNA-Seq analysis. After stringtie using ballgown I get FPKM and TPM values for every gene.
I have seen that edgeR, Deseq2 can be used for Counts data. I would like to know which R package needs to be used for differential analysis with TPM values?
Which one is better for differential analysis FPKM or TPM?