Question: Extracting genes from one file using coordinate from second file
0
gravatar for bhuwankhatri11
8 days ago by
bhuwankhatri1110 wrote:

Hi Alex Reynolds,

I am sorry to ask this question once again. Because I got some error from previous command line.

My first file looks like: First.txt

deletion    chr10:1726501-1755000   28500   0.586226    9.73037E-05 715.754 0.00171548  3546.87 0.114216    1.17241

Second file looks like: Second.txt

Chr10   NC_029525.1 gene    1672245 1676954 -   LOC107318572
Chr10   NC_029525.1 gene    1677076 1682931 -   C10H15orf39
Chr10   NC_029525.1 gene    1690899 1710413 -   PPCDC
Chr10   NC_029525.1 gene    1710723 1714472 -   LOC107318577
Chr10   NC_029525.1 gene    1714558 1714977 -   LOC107318579
Chr10   NC_029525.1 gene    1717116 1719122 +   RPP25
Chr10   NC_029525.1 gene    1721742 1725395 +   LOC107318578
Chr10   NC_029525.1 gene    1725935 1728167 +   FAM219B
Chr10   NC_029525.1 gene    1728336 1731151 -   MPI
Chr10   NC_029525.1 gene    1731194 1739576 +   LOC107318570
Chr10   NC_029525.1 gene    1739821 1743801 +   ULK3
Chr10   NC_029525.1 gene    1744568 1747749 -   CPLX3
Chr10   NC_029525.1 gene    1752411 1759515 -   CSK

I want to get result like this: Asnwer.txt

Chr10   NC_029525.1 gene    1725935 1728167 +   FAM219B
Chr10   NC_029525.1 gene    1728336 1731151 -   MPI
Chr10   NC_029525.1 gene    1731194 1739576 +   LOC107318570
Chr10   NC_029525.1 gene    1739821 1743801 +   ULK3
Chr10   NC_029525.1 gene    1744568 1747749 -   CPLX3
Chr10   NC_029525.1 gene    1752411 1759515 -   CSK
bedoops awk bash • 152 views
ADD COMMENTlink modified 8 days ago by Alex Reynolds21k • written 8 days ago by bhuwankhatri1110
1
gravatar for Alex Reynolds
8 days ago by
Alex Reynolds21k
Seattle, WA USA
Alex Reynolds21k wrote:

Preprocess the first file into a sorted BED file.

You can split the second field of each line by a regular expression pattern and write the output to a three-column, unsorted BED file, sort it, and write the sorted output to a BED file:

$ awk -v OFS="\t" '{ n=split($2,a,/[:-]/); print a[1], a[2], a[3]; }' First.txt | sort-bed - > First.bed

Preprocess the second file into a sorted BED file.

This step requires decapitalizing the chromosome name (so that you can do set operations later), shuffling columns into UCSC BED order, and sorting:

$ awk -v OFS="\t" '{ if ($1~/^Chr/) { $1=tolower(substr($1,1,1))substr($1,2); } print $1, $4, $5, $2, $3, $6, $7; }' Second.txt | sort-bed - > Second.bed

Double-check that you have BED files by running head First.bed and head Second.bed. These should look like sorted BED files. If not, review the steps before proceeding.

Now you can do set operations, filtering the Second.bed file based on overlaps with elements in First.bed:

$ bedops --element-of 1 Second.bed First.bed > Filtered.bed

Re-capitalize the chromosome name in Filtered.bed and re-order columns to match your desired format:

$ awk -v OFS="\t" '{ $1=toupper(substr($1,1,1))substr($1,2); print $1, . . . }' Filtered.bed > Filtered.txt

I'll let you fill in . . . so that you get a little experience with moving columns around with awk.

ADD COMMENTlink modified 8 days ago • written 8 days ago by Alex Reynolds21k

Thank you Alex for quick reply. I should say it worked partially because it is missing some genes list from chromosomes like chrZ, chrM, LGE64 and LGE22C19W28_E50C23.

ADD REPLYlink written 8 days ago by bhuwankhatri1110

Yep, that's a bug. See the revised preprocessing step for the second file.

ADD REPLYlink written 8 days ago by Alex Reynolds21k

I am able to get chrZ but chromosomes chrW, chrM, LGE64 and LGE22C19W28_E50C23 are still missing.

ADD REPLYlink modified 8 days ago • written 8 days ago by bhuwankhatri1110

I added a fix to test for the chromosome name in the second preprocess step. If the name starts with Chr then the first letter is made lowercase. Otherwise, the name is left unmodified.

You'll need to investigate further on your own, from this point out, as you may need to preprocess your first file, as well, depending on how you have named chromosomes in that file.

ADD REPLYlink written 8 days ago by Alex Reynolds21k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1353 users visited in the last hour