grep sorting question
Hi I have the following code:
grep -w -F -i -f gene_list.txt gencode.vM18.promoters.bed > gene_list_promoters.bed
head gene_list.txt
Rn45s
Malat1
Sptbn1
head gencode.vM18.promoters.bed
chr1 3071252 3075252 ENSMUSG00000102693.1 . + HAVANA gene . ID=ENSMUSG00000102693.1;gene_id=ENSMUSG00000102693.1;gene_type=TEC;gene_name=RP23-271O17.1;level=2;havana_gene=OTTMUSG00000049935.1
chr1 3100015 3104015 ENSMUSG00000064842.1 . + ENSEMBL gene . ID=ENSMUSG00000064842.1;gene_id=ENSMUSG00000064842.1;gene_type=snRNA;gene_name=Gm26206;level=3
chr1 3669498 3673498 ENSMUSG00000051951.5 . - HAVANA gene . ID=ENSMUSG00000051951.5;gene_id=ENSMUSG00000051951.5;gene_type=protein_coding;gene_name=Xkr4;level=2;havana_gene=OTTMUSG00000026353.2
It retrieves the promoters specific to the .txt list of genes, however this is an ordered gene list and the promoters that are retrieved are ordered by chromosome, wondering what I can add to the code the order the promoters in the same order as gene list? Thanks.
1 answer
What is the purpose of this exercise?
I think the best way to do this would be to use R, especially if the number of lines in the two files are not the same. If you're willing to abuse bash and re-create the file, you can loop through gene_list.txt line by line and run the grep for each line, >>-ing to the output file. But like I said, that would be abusing the shell.
Log in to answer this question.
More posts like this
-
featurecounts not working on mirbase annotation file
written by demoraesdiogo2017 •
-
Is there a way to safely remove all underscore from gene id for hg19 from genecode
written by simplitia •
-
Error with Rseqc infer_experiment - unknown datatype
written by newbie
-
I have a list of chromosomes, positions, gene ids, and transcript ids (from a gencode GTF file) - h…
written by claudiadast •
-
How to safely rename the chromosome names in fasta and gtf
written by ddzhangzz
-
Choice of Gene Length in TPM Calculations
written by dec986
-
extracting a specific field from .bed file and sort into new file
written by g.sathish.k •
-
Comprehensive gene annotation for drosophila melanogaster
written by newscient •
-
Cufflinks seg fault
written by shg018 •
-
extract only geneID and gene symbol from GTF file
written by catherine
You should add the output of
headfor each of the two input files.In all probability, you'll either need to do the sorting in a separate step or use process substitution.
Just added the
headoutputWould this work?
Update:
instead of
grep ${line}please usegrep "gene_name=${line};"or something similar.Also If you are looking to extract based on gene_list.txt I think your code has a bug, because I'm only testing on three first genes and first 3 lines of gencode.vM18.promoters.bed and even though these genes are not found in the toy bed file your code outputs something.
Toy Gene list
Toy bed file
Your code:
Your grep will pick the entry for
A1BG-AS1when grepping forA1BG.Thanks for pointing this out.
I think
grep "gene_name=${line};"should fix this.This does seem to work in terms of ordering the sites however the gencode list has 54K rows and my gene list has 14k but I get a final tally of 140K rows. Is that due to alternate promoters?
>>appends, so each time you run the command it adds to the same file. You might want to delete the resultfile (gene_list_promoters.bed ) and run again because the result of the first run is not correct.Also please make sure your gene_list.txt doesn't have duplicate names. Use
To make sure it's equal the number of lines in gene_list.txt