Tuesday, November 4, 2014

GC contents with BEDtools

You can calculate the GC content of your sequence with BEDTools. (source)

Usage: bedtools nuc [OPTIONS] -fi -bed [FILENAME]


For example, If you have your sequence file named "file.fasta" as a FASTA format, first you need to generate a text file named "interval.bed" with position information of your interests. (Of course your interval can be whole genome or chromosome)

  chr1 1000 2000
  chr2 1000 2000
  ....

Then execute following command:


  bedtools nuc -fi file.fasta -bed interval.bed


Output might look like below:

#1_usercol 2_usercol 3_usercol 4_pct_at 5_pct_gc 6_num_A 7_num_C 8_num_G 9_num_T 10_num_N 11_num_oth 12_seq_len
      chr1     1000     2000 0.560000 0.440000 355 253 187 205 0 0 1000


First three columns (usercol) show7 your BED file entry, followed by information generated by your bedtool command as below.

Output format:
  The following information will be reported after each BED entry:
    1) %AT content
    2) %GC content
    3) Number of As observed
    4) Number of Cs observed
    5) Number of Gs observed
    6) Number of Ts observed
    7) Number of Ns observed
    8) Number of other bases observed
    9) The length of the explored sequence/interval.
    10) The seq. extracted from the FASTA file. (opt., if -seq is used)
    11) The number of times a user's pattern was observed.
        (opt., if -pattern is used.)


You can also print sequence of DNA in intervals as below:


  bedtools nuc -seq -fi file.fasta -bed interval.bed


and can even look for specific pattern in that sequence:


  bedtools nuc -seq -pattern 'GATACA' -fi file.fasta -bed interval.bed


No comments:

Post a Comment