Wednesday, January 30, 2013

How to use grep with files located within subfolders

If you want to search multiple files for lines with a specific string within all the subfolders of a current folder use grep with -r option:

grep -r "string" .

Wednesday, January 23, 2013

Concatenate files in Unix

To concatenate two files in Unix use cat:
cat file1.txt file2.txt > output.txt

Copy specific files from all subfolders of a folder to another folder in Unix

If you want to copy all the .txt files for example from the subfolders of the current folder to another destination folder do the following:

find . -name *.txt -type f -exec cp -v {} ./all_txt_files/ \;

This will copy all the txt files in a current folder and in all subfolders of a current folder to the folder all_txt_files located in the current folder.

However, note that this will only look in the subfolders and not sub-subfolders and so on.

To look in the sub-subfolders as well simplyput the name of the file you are searching in parenthesis:

find . -name "*.txt" -type f -exec cp -v {} ./all_txt_files/ \; 

And also if you have the files with the same name they will overwrite each other,  so only the last file copied remains.

Tuesday, January 22, 2013

How to select specific genes from whole gene lists using grep in Unix

If you have a big list of genes and want to select specific ones you can easily do it using the grep command, e.g.

grep "Fgf5" gene_exp.diff
will select for all the lines containg Fgf5 gene name

Now if you have a gene that belongs to a gene family, there may be a possibility where you will using this command pick up the genes you don't want. e.g.

grep "Adamts1" gene_exp.diff
will give as output Adamts1, but also Adamts10, Adamts11, Adamts12 lines etc.
In this case you may use this specific trick. If the next column in file contains the coordinates of the gene (e.g. chr17:33661049-33690727), as it is in the gene_exp.diff file made by Cuffdiff, do the following:

grep "Adamts1.c" gene_exp.diff
Here .c represents that one character is separating Adamts1 from the c character from the gene coordinates, which is only the case for Adamts1 lines, and not for Adamts10 lines etc. where there are two characters until c charachter.

Monday, January 21, 2013

Remove without confirmation using rm command in Unix

If you have multiple files that you want to remove using rm command in Unix, e.g.
rm *.txt
You will be prompted to confirm each of the files.
This may take long time if you have a lot of files, so IF YOU ARE CERTAIN YOU WANT TO DELETE THESE FILES (SINCE THIS CANNOT BE UNDONE) do the following:

rm -frv

-f Attempt to remove the files without prompting for confirmation, regardless of the file's permissions.
-v Be verbose when deleting files, showing them as they are removed.
-r Equivalent to -R. Recursive. To include subfolders.

How to make custum labels for the axis in R

If you want to make your own labels for the x- or y- axisdo the following:
Delete the x- or y-axis with xaxt='n', or yaxt='n'.
Then make a vector containing the labels for your axis, e.g.
fold.change <- c("4","3","2","1.5","1.2","all")
Then make the axis:
axis(side = 1, at=1:6, labels = fold.change, tck = -0.02)
side=1 will define the side of the axis, at=1:6 will split the axis into 6 parts, labels=fold.change will put the labels from the vector fold.change, tck=-0.02 will define the orientation and size of the ticks.


The following code will produce the graph bellow.

> plot(x, main="Differentially expressed genes, RN vs KD, p<0.05", xlab="Fold change", ylab="Number of genes", col="red", xaxt='n', type ="b")
> points(y, col="blue", type="b")
> points(z, col="black", type="b")
> points(a, col="orange", type="b")
> fold.change <- c("4","3","2","1.5","1.2","all")
> axis(side = 1, at=1:6, labels = fold.change, tck = -0.02)
> mtext("Upregulated - red, downregulated - black. Filtered with RPKM>1 in blue and orange")






How to plot two or more variables on the same graph in R

If you want to plot two or more variables in the same graph in R, use the following code:
First plot one variable using the standard plot command.
The second variable and every next one you can place on the graph using the points command.
The graph bellow was made using the following code:

> plot(x, main="Upregulated genes, RN vs KD, p<0.05", xlab="Fold change", ylab="Number of genes", col="red", xaxt='n', type ="b")
> points(y, col="blue", type="b")
> fold.change <- c("4","3","2","1.5","1.2","all")

> axis(side = 1, at=1:6, labels = fold.change, tck = -0.02)
> mtext("Data with filtering RPKM>1 in blue")









Beware not to overwrite output files of Cuffdiff by mistake

I have set Cuffdiff to process several files with the following command:

C02HQ105DHJQ:Lluis - EB project mpjanic$ cuffdiff merged_asm/merged.gtf -o diff_out -b genome.fa -p 2 --no-diff --labels 4_BMI1_3604_ACAGTG,4_DCAF7_3605_GTGAAA,4_EED_3610_ACTGAT,4_MEL18_3602_GCCAAT,4_PCGF5_3603_CTTGTA,4_RING1B_3606_GTCCGC,5_CBX7_3612_CTTGTA,5_CTR_3611_GCCAAT,ES-CBX7_3559_CTTGTA,ES-CTR_3558_GCCAAT,ES-MEL18_3561_GTCCGC -u \
> ./4_BMI1_3604_ACAGTG/accepted_hits.bam \
> ./4_DCAF7_3605_GTGAAA/accepted_hits.bam \
> ./4_EED_3610_ACTGAT/accepted_hits.bam  \
> ./4_MEL18_3602_GCCAAT/accepted_hits.bam \
> ./4_PCGF5_3603_CTTGTA/accepted_hits.bam \
> ./4_RING1B_3606_GTCCGC/accepted_hits.bam \
> ./5_CBX7_3612_CTTGTA/accepted_hits.bam \
> ./5_CTR_3611_GCCAAT/accepted_hits.bam \
> ./ES-CBX7_3559_CTTGTA/accepted_hits.bam \
> ./ES-CTR_3558_GCCAAT/accepted_hits.bam \
> ./ES-MEL18_3561_GTCCGC/accepted_hits.bam


After the extensive run has finished I have noticed the output files were empty:

C02HQ105DHJQ:diff_out mpjanic$ ls -li
total 0
2605038 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 cds.count_tracking
2605031 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 cds.diff
2605034 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 cds.fpkm_tracking
2605042 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 cds.read_group_tracking
2605028 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 cds_exp.diff
2605027 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 gene_exp.diff
2605039 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 genes.count_tracking
2605035 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 genes.fpkm_tracking
2605043 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 genes.read_group_tracking
2605025 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 isoform_exp.diff
2605036 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 isoforms.count_tracking
2605032 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 isoforms.fpkm_tracking
2605040 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 isoforms.read_group_tracking
2605030 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 promoters.diff
2605044 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 read_groups.info
2605045 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 run.info
2605029 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 splicing.diff
2605026 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 tss_group_exp.diff
2605037 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 tss_groups.count_tracking
2605033 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 tss_groups.fpkm_tracking
2605041 -rw-r--r--  1 mpjanic  staff  0 Jan 21 03:07 tss_groups.read_group_tracking



I have noticed that after the run was finished by mistake this command was executed (most probably by me typing the Up arrow and [, followed by Enter):

C02HQ105DHJQ:Lluis - EB project mpjanic$ cuffdiff merged_asm/merged.gtf -o diff_out -b genome.fa -p 2 --no-diff --labels 4_BMI1_3604_ACAGTG,4_DCAF7_3605_GTGAAA,4_EED_3610_ACTGAT,4_MEL18_3602_GCCAAT,4_PCGF5_3603_CTTGTA,4_RING1B_3606_GTCCGC,5_CBX7_3612_CTTGTA,5_CTR_3611_GCCAAT,ES-CBX7_3559_CTTGTA,ES-CTR_3558_GCCAAT,ES-MEL18_3561_GTCCGC -u ./4_BMI1_3604_ACAGTG/accepted_hits.bam ./4_DCAF7_3605_GTGAAA/accepted_hits.bam ./4_EED_3610_ACTGAT/accepted_hits.bam  ./4_MEL18_3602_GCCAAT/accepted_hits.bam ./4_PCGF5_3603_CTTGTA/accepted_hits.bam ./4_RING1B_3606_GTCCGC/accepted_hits.bam ./5_CBX7_3612_CTTGTA/accepted_hits.bam ./5_CTR_3611_GCCAAT/accepted_hits.bam ./ES-CBX7_3559_CTTGTA/accepted_hits.bam ./ES-CTR_3558_GCCAAT/accepted_hits.bam ./ES-MEL18_3561_GTCCGC/accepted_hits.bam[
You are using Cufflinks v2.0.2, which is the most recent release.
open: No such file or directory
File ./ES-MEL18_3561_GTCCGC/accepted_hits.bam[ doesn't appear to be a valid BAM file, trying SAM...
Error: cannot open alignment file ./ES-MEL18_3561_GTCCGC/accepted_hits.bam[ for reading


This has overwritten the output files. I have checked this by modifing gene_exp.diff file to contain 5 bytes and then running this comand again and after this the file contained again 0 bytes.

So be careful not to touch the arrows when executing Cuffdiff or any similar tool in Unix.