How to generate random genomic sequences that are GC- and length-matched to my input sequences?

Question

[Xianjun]

Hi,

I want to generate a list of random genomic sequences that match GC-content and length to my input DNA sequences. Note that my input sequences are not of the same length. My expected output will be the same number of random sequences that match length and GC-content to the individual input sequence.

>seq1
AACAATGAGTTGTCATTTTATCCAAATCTAAAAAAAAACATACATACAGTATCAAGTTCTGGTTAAAGTAGTAGATGTACTTTACGTGGTACATAAGAGTATCTATACAGTTTGTGGG
>seq2
atacaacctgaaaaacagagagggaaaaataattaagaaaaatgaatagaagcttgggacatccttaatgagagtaccagaggagatgtgagagagaaagaagcagaaaaaaagttccaagcagaaaaaatgtccaagaaataatagctgaaaagctcccaaatttgctgaaaaatgttaacctacacattcaagaagcccaaaaagctctacacaaagagacacacctagacacacacacctggacacacacacctgga
>seq3
cacctaagctggagtgcagtggtgcaattacagcttactgcagcctcaatctcccaggcttaagggatcctcccatgtagctgagactacaggcatgagccactatgtccagctaatttttaaattttttgtagagacagggtctcgctaccttgaacgggctgaccttgaactcctgggctctggtggccctcctgtgttgacctcccaaagcattgggattacaggcatgagccactgcacccagccTAGAAGCTCTGTTCATATTTATTT
>seq4
CATGAGGCCCAGTCTGTGAACAGAGACCAGGTCTAACCCCTTCTTCCAGGAAAGCCTCGTAGGGCCTTCTGGCCAAGAGGCCACGagtggtgaagactgcagactctgaaatcagaaatacctgggctccactgtcagcatggcagctgaggaagagtgaaaattcctctaagttcttttagaagtcccagcctccccatgtaactggggaACTGATGGGAGGAGCAGAGCTGTCTGTGCACATAAGAAGTTCTCAGTAAATGGAGACAGTTACTATTTCTGTTATTATTGAATTTGAACAAATTCCCTGGGTATGTGTGGGGGGACACTTCAGGTGAAAACACGCCCCTCCTCCCCTGGTGCGGGGGCCTGTGCTGCCACCCTCTGGAAGCCTGCAGAGGGGCAGGGAAAACAGACCCTGAACAAAAGTGTGCACCCAGTGAGGAGGTGCAAGGGCACAAAGGTGGCACCAAGTGCCTCAAGGAGAGGCTGAAACGCGGCCTGGGGACCTCGCAGTGGTCTGGTCATATAGGC
>seq5
GGCATGTGGTGTCAGCAGAGGTGCCTCAAGGATAGAGTGAGTCCAGAGTCTAGAAAGGAGCAGATCACCAGGCTCTGGGAAGAGCACAGCATGGGTGCACACACTGCTCTACCCAGCATGGCTGCCGACCCAAAGACAGCAAAGCCAAGAAGGACACACAAGCGTGGCCAGATGCAGCCCTGTGAGGAAACTTACCCAAGAACGGGACGATGGGCTTGAGAAACCatccatctacaaggatggcgtttgctgcagcaatgtttataataaattgtgggaaactgtgaactgcctaaatgtctcacaataggaacaaattagtgcaccacaccatgaaactctctacagcTCCTGAGTTACAGAACGACAGTATAATACTA

Could anyone please suggest a proper tool for this purpose?

btw, I have looked into NullSeq, GC_compo, and BiasAway. But none of them fits my needs.

Thanks,

-Xianjun

[GenoMax]

Take a look at randomreads.sh from BBMap suite. It should produce fasta files.

[Xianjun]

Thanks for the reply. I looked into the option of randomreads.sh, right it can produce fasta files, but it doesn't match the length distribution and GC content.

[GenoMax]

You may need to write something yourself perhaps since you have a specific need.

[Xianjun]

Sure. I thought there might already be some tools available since it's a common task when people are asked to compare to a length- and GC-matched background set, e.g. TFBS enrichment analysis.

Answer 1

[finswimmer]

The easiest way could be to read in every sequence and shuffle the characters.

fin swimmer

[GenoMax]

Simple solutions are the best. Depending on input from OP you can move this to an answer.

[Xianjun]

yes, that's also what I have in my mind so far. I am wondering if any "better" background (e.g. part of genomic region). I like the idea of http://opossum.cisreg.ca/GC_compo/ in oPOSSUM, where they use background regions as human open chromatin regions compiled from the lowest scoring ChIP-Seq peaks in the Encode project. The short thing is their output sequences have a fixed length, not matching my input length distribution.

[WouterDeCoster]

I did something similar for fastq format, should give some ideas on how to adapt this to fasta: https://github.com/wdecoster/fastq-anonymous

Answer 2

[dariober]

Among the many programs that I hacked together, there is MarkovSequenceGenerator which may be you can make use of. As I remember, just like the OP, I needed a null distribution of sequences matching nucleotides and di-nucleotide content. So this program shuffles sequences in a fasta file using a Markov model so that the shuffled sequences look like the original one.

wget https://github.com/dariober/Java-cafe/releases/download/v0.1.0/MarkovSequenceGenerator-0.1.0.jar

java -jar MarkovSequenceGenerator-0.1.0.jar -h

usage: MarkovChainGenerator [-h] -f FASTA [-nr NRANDOM] [-l LENGTH] [-o OUTPUT] [-n NORDER]
                            [-p PCT_PSEUDOCOUNTS] [--version]

DESCRIPTION
Generate random sequences from input sequence using Markov chain model

optional arguments:
  -h, --help             show this help message and exit
  -f FASTA, --fasta FASTA
                         Input fasta file to generate sequences from. Can be gzipped.
  -nr NRANDOM, --nrandom NRANDOM
                         Number of chains to generate. *Each* sequence in fastainput will be randomized this many times (default: 1)
  -l LENGTH, --length LENGTH
                         Lenght of each random sequence. Default to length of input sequence (default: -1)
  -o OUTPUT, --output OUTPUT
                         Output file. Use - to stdout (default: -)
  -n NORDER, --norder NORDER
                         Order of the Markov model. Must be greater than 0. (default: 1)
  -p PCT_PSEUDOCOUNTS, --pct_pseudocounts PCT_PSEUDOCOUNTS
                         If the input sequence does not contain all the possible  kmers,add  the  missing  ones  at  this percentage of the total counts. E.g. -p 0.1
                         makes all the missing kmers to sum up to no more  than  10%  of  the  total  actual counts. Setting to zero might terminatethe chain with an
                         error. (default: 0.1)

Answer 3

[Hughie]

I wrote a pure bash script for this purpose, verbose but useful. Wish to see better versions :)

Negetive_sequence_matched_length_GC.sh

This script is for generating length and GC content matched background sequence, heavily using bedtools
Negetive_sequence_matched_length_GC.sh [input bed file] [genome size] [genome fasta] [iteration times]
[input bed file]: A bed file as template, just first 3 cols are used
[genome size]: Genome size for bedtools
[genome fatsa]: Genome Fatsa file for bedtools
[iteration times]: Sequence pools for comparison, 100 is pretty good