MY SCRIPT (11): GC analysis of genes ........

#To generate the data file
#! usr/bin/python
#Execute as: GC_content.py
from Bio import SeqIO
#To find the count of A, T, C, G  9their % too) in fasta sequence
#Opened a FASTA file
input_file = open('/home/pseema/denovo_analysis/genome_ffn_files/L_12152015.ffn', 'r')
#Opened a output file
output_file = open('GC_content','w')
#Header line
#Header looks likie Gene    A    C    G    T    Length    CG%
output_file.write('Gene\tA\tC\tG\tT\tLength\tCG%\n')


for cur_record in SeqIO.parse(input_file, "fasta") :
#count nucleotides in this record...
    gene_name = cur_record.name
    A_count = cur_record.seq.count('A')
    C_count = cur_record.seq.count('C')
    G_count = cur_record.seq.count('G')
    T_count = cur_record.seq.count('T')
    length = len(cur_record.seq)
    cg_percentage = float(C_count + G_count) / length
    output_line = '%s\t%i\t%i\t%i\t%i\t%i\t%f\n' % \
    (gene_name, A_count, C_count, G_count, T_count, length, cg_percentage)
    output_file.write(output_line)

output_file.close()
input_file.close()

#To manipulate the data file
#! usr/bin/bash
#run as:  sh GC_analysis.sh
#To process GC_content datafile generated bt GC_content.py
#Run as python GC_content.py
#mkdir /home/pseema/denovo_analysis/result_files/GC_data

#There are 4423 lines. Remove lines starting with 'gnl'. It has 4173 lines
 sed '/^gnl/d' GC_content > /home/pseema/denovo_analysis/result_files/GC_data/truncated_GC_content

 #Extract headers of the original .ffn file
 grep "^>" /home/pseema/denovo_analysis/genome_ffn_files/L_12152015.ffn > /home/pseema/denovo_analysis/result_files/GC_data/L_12152015.only_header |head

 #There are 4422 lines. Remove lines starting with 'gnl'. It has 4172 lines
 sed '/^>gnl/d' /home/pseema/denovo_analysis/result_files/GC_data/L_12152015.only_header  > /home/pseema/denovo_analysis/result_files/GC_data/truncated_L_12152015.only_header

#Remove  line 1 of truncated_GC_content
awk 'NR!=1' /home/pseema/denovo_analysis/result_files/GC_data/truncated_GC_content > /home/pseema/denovo_analysis/result_files/GC_data/no_legend_GC_content

#Take all columns of this file and add side by side to truncated_L_12152015.only_header file
#Extract all columns except 1
awk '{$1=""; print $0}' /home/pseema/denovo_analysis/result_files/GC_data/no_legend_GC_content > /home/pseema/denovo_analysis/result_files/GC_data/all_fields_except_first
#Paste side by side with truncated_L_12152015.only_header
paste /home/pseema/denovo_analysis/result_files/GC_data/truncated_L_12152015.only_header  /home/pseema/denovo_analysis/result_files/GC_data/all_fields_except_first > /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields

#Add a legend and format.
{ printf 'Gene\tProtein_name\tA\tC\tG\tT\tLength\tGC_percentage\n'; cat /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields;} > /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend
#Format.Formatting is not working because the protein name is spanning multiple words
#column -t /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/formatted_pasted_fields.tsv



grep " PE family " /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/PE_only
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/PE_only > /home/pseema/denovo_analysis/result_files/GC_data/PE_only_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/PE_only_percentage > /home/pseema/denovo_analysis/result_files/GC_data/PE_only_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/PE_only_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/PE_only_percentage_sorted
echo "No. of PE genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/PE_only_percentage_sorted

grep " PPE family " /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/PPE_only
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/PPE_only > /home/pseema/denovo_analysis/result_files/GC_data/PPE_only_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/PPE_only_percentage > /home/pseema/denovo_analysis/result_files/GC_data/PPE_only_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/PPE_only_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/PPE_only_percentage_sorted
echo "No. of PPE genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/PPE_only_percentage_sorted


grep " hypothetical " /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/hypothetical
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/hypothetical > /home/pseema/denovo_analysis/result_files/GC_data/hypothetical_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/hypothetical_percentage > /home/pseema/denovo_analysis/result_files/GC_data/hypothetical_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/hypothetical_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/hypothetical_percentage_sorted
echo "No. of hypothetical genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/hypothetical_percentage_sorted



grep "phage\|Phage" /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/phage
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/phage > /home/pseema/denovo_analysis/result_files/GC_data/phage_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/phage_percentage > /home/pseema/denovo_analysis/result_files/GC_data/phage_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/phage_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/phage_percentage_sorted
echo "No. of phage genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/phage_percentage_sorted


grep "transposase\|Transposase" /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/transposase
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/transposase > /home/pseema/denovo_analysis/result_files/GC_data/transposase_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/transposase_percentage > /home/pseema/denovo_analysis/result_files/GC_data/transposase_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/transposase_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/transposase_percentage_sorted
echo "No. of transposase genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/transposase_percentage_sorted



grep "toxin\|Toxin" /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/toxin
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/toxin > /home/pseema/denovo_analysis/result_files/GC_data/toxin_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/toxin_percentage > /home/pseema/denovo_analysis/result_files/GC_data/toxin_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/toxin_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/toxin_percentage_sorted
echo "No. of toxin genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/toxin_percentage_sorted



grep "methyltransferase\|Methyltransferase" /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase > /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase_percentage > /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase_percentage_sorted
echo "No. of methyltransferase genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/methyltransferase_percentage_sorted

grep ">L_0" /home/pseema/denovo_analysis/result_files/GC_data/pasted_fields_with_legend > /home/pseema/denovo_analysis/result_files/GC_data/all
awk '{print $1, $2, $3, $4, $NF}' /home/pseema/denovo_analysis/result_files/GC_data/all > /home/pseema/denovo_analysis/result_files/GC_data/all_percentage
#Multiply column 5 with a value (value is 100 here)
awk 'BEGIN{FS=OFS=" "}{print $1,$2,$3,$4,$5*a}' a=100 /home/pseema/denovo_analysis/result_files/GC_data/all_percentage > /home/pseema/denovo_analysis/result_files/GC_data/all_percentage_multiplied
cat  /home/pseema/denovo_analysis/result_files/GC_data/all_percentage_multiplied| sort -k5,5 > /home/pseema/denovo_analysis/result_files/GC_data/all_percentage_sorted
echo "No. of all genes......"
wc -l /home/pseema/denovo_analysis/result_files/GC_data/all_percentage_sorted