Assigning gene to a chromosome segment from CBS segmentation output:

[angerusso]

Hi All,

I have an CBS segmentation algorithm output for 10 tumor samples each from 2 different tumors.

Now, I am in an urgent need to assign gene (followed by all genes) that belong to a particular chromosome segment after I removed all the CNVs from segment data. The format of the data is:

Sample Chromosome Start End Num_Probes Segment_Mean
Sample1A-TA 1 51598 76187 15 -1.115

Could anyone suggest an R library or code or method that I can quickly use to get the genes assigned to CBS output.

Thanks so much,
Angel

[zee]

Hi

You could convert your segment data to bed first and then match genes with Bedtools.
Start by converting your data to bed format:

perl -lane 'printf "chr$F[1]\t%s\t$F[3]\t$F[0]\t$F[5]\n",$F[2]-1'  segments.txt > segments.bed

Open a new browser window here and click "Get output"

Install bedtools.

Unzip the bed file and match your regions to the genes:

intersectBed -wa -wb -a  segments.bed -b  refseq.bed > result.txt

You should have a line of matching refseq for regions matching these in the genome

[angerusso]

Hi Zee,

Thanks very much for your reply and I used BedTools as you suggested. One quick question:

I am supposed to report the format of the file as:

sample chr start end genestart geneend #genesinsegment namesofgenesinsegment

Now I need to know how "intersectbed" is finding the intersection. I also saw there are options to give the amount of overlap using "-f" option. Should I play with these different amounts of overlap. Are we assuming that all genes that are within the boundaries of the chromosome are the only ones reported?

Output from one segment leads to 2 genes as below

The output is from above segment is two genes in the segment:
chr1 51598 76187 Sample1A-TA 15 -1.115 chr1 52878 53750 ENSG00000205292 AL627309.15-4 pseudogene
chr1 51598 76187 Sample1A-TA 15 -1.115 chr1 58954 59871 ENSG00000177693 OR4F5 protein_coding

I can just simply write a perl code to get the format I want and I still will be reporting two lines per gene, as there are two genes in that segment.

Cheers,
Angel

[zee]

Hi Angel,

The Bedtools PDF manual has some very clear examples and illustrations of how the different functions work, including intersectBed. Also have a look at the filo package for summing up your data. There is a groupBy function in filo that will collapse the list of genes for you.

Use the output of BedTools to make the format you need.

[angerusso]

Hey,
Thanks so much!

Angel