Converting Chip-seq fastq format ({$x_k$} list of seq segments) to genome browser {sequence $x_i$, count-number $y_i$} format

Question

Very new to bioinformatics. I am trying to figure out the steps required to convert a chip-seq fastq file into a genome browser track, with the end goal of being able to access the sequence (x) vs. activity level/count (y) data associated.

Any genome browser activity/count track data has the following functional format $y=F(x)$

Here x is the sequence axis, and y is the activity/count axis. Typically the x(i) is the nucleotide associated with the precise genomic coordinate i (an integer number) depending on the genome used to align (human, mouse, etc.). What I need is a format allowing direct access to the nucleotide representation x(i), associated with the index i, as well as the count levels at y(x(i)). What is the closest data file type that I need? And how do I get to it?

In other words, what are the steps for a fastq file associated with a chip-seq data to be transformed to this x(i),y(i) format, or best yet if there is an example Python or R code out there that could illustrate the steps of going from file1.fastq to something like to two files file1_x_sequences and file2_y_counts? I only know the first step is an alignment to a genome, but stuck on the next steps to complete the picture. Thank you.

Ultimately I need a list of sequences (tags) lets say 250bp long, and a count number that represents the chip-seq enrichment for the part of the genome (the 250bp region). Further reading, indicates I need to detect/call "peaks" with MACS2/homer to generate a BAM file and ultimately a Fasta file for the sequence lists, but I don't know yet how the 'counts' come in (in what file format and where from). Any hints there would complete the picture I think. Ultimately this is training data for x=250bp sequence to y=count/chip-seq enrichment activity model.

Answer 1

I have no idea about what you are actually trying to achieve, but getting these browser tracks is simple. First you align your fastq files to a reference genome, e.g. with bowtie2, this returns a SAM file. Sort the SAM file, convert to BAM and sort + index it with samtools. Then use something like bamCoverage from the deeptools suite to get your browser tracks.

A word of caution. Since you seem to be new to OMICS and bioinformatics in general: NGS data, especially applications that are as noisy as ChIP-seq have lots of pitfalls. Be sure to go through the literature to find references for what you are trying to develop, at least the basic idea and discuss with experienced folks whether it makes sense. It sounds like you want counts of the peak summits, is that correct? In any case, be sure not to reinvent the wheel and to do a meaningful analysis. Discussion with experienced peers is key, that may save you from pitfalls and rookie mistakes. This entire choice of words "sequence axis", "activity level" is very uncommon, I have never heard of this before despite working on NGS data for years, for what it's worth.