FASTQ HANDS-ON TUTORIAL 
====================================================
To download dataset
----------------------------
wget https://ftp.sra.ebi.ac.uk/vol1/fastq/ERR105/000/ERR10549500/ERR10549500_1.fastq.gz

----------------------------------------------------
1. calculate read-length distribution
----------------------------------------------------

zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {print length($0)}' | \
sort | uniq -c

---------------------------------------
AVG read length
------------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {sum += length($0); n++} END {print sum/n}'

---------------------------------------
MIN read length
------------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {print length($0)}' | \
sort -n | head -1 

---------------------------------------
MAX read length
------------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {print length($0)}' | \
sort -n | tail -1


----------------------------------------------------
GC content
----------------------------------------------------
----------------------------
For first few sequences
------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {print $0}' | head

--------------------------------------------
Calculate GC percentage for each read :
-----------------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {
gc=gsub(/[GCgc]/,"");
total=length($0);
print (gc/total)*100
}' | head

--------------------------------
AVG GC percentage 
--------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {
gc=gsub(/[GCgc]/,"");
total=length($0);
sum += (gc/total)*100;
n++
}
END {
print "Average GC% =",sum/n
}' 


--------------------------------
Count A,T,G,and C Nucleotide composition
----------------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2' | \
tr -d '\n' | \
fold -w1 | \
sort | uniq -c

-----------------------------
Count #nucleotides per read for all reads
--------------------------------------------
zcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {
    read++;
    A=gsub(/A/,"");
    T=gsub(/T/,"");
    G=gsub(/G/,"");
    C=gsub(/C/,"");
    N=gsub(/N/,"");

    printf "Read_%d: A=%d T=%d G=%d C=%d N=%d\n",
           read,A,T,G,C,N
}'

---------------------------------------
Count #nucleotides per read first 10 reads
-------------------------------------------
gzcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {
    read++;

    A=gsub(/A/,"");
    T=gsub(/T/,"");
    G=gsub(/G/,"");
    C=gsub(/C/,"");
    N=gsub(/N/,"");

    printf "Read_%d: A=%d T=%d G=%d C=%d N=%d\n",
           read,A,T,G,C,N;

    if (read == 10) exit;
}'


------------------------------------------------
Search for N bases : Count reads containing N: A ?

-------------------------------
gzcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 && /N/ {count++} END {print count}' 

---------------------------------------------
Percentage of reads that contain at least one N: 
----------------------------------------------------
gzcat ERR10549500_1.fastq.gz | \
awk '
NR%4==2 {
total++;
if ($0 ~ /N/) n++;
}
END {
print "Reads with N =",n;
print "Total reads =",total;
print "Percentage =",100*n/total,"%"
}' 

----------------------------
To generate read_length.txt
-----------------------------
gzcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {print length($0)}' \
> read_lengths.txt 


----------------------------------
Python code for plot_histogram.py
-----------------------------------
import matplotlib.pyplot as plt
# Read the read lengths
with open("read_lengths.txt") as f:
    lengths = [int(line.strip()) for line in f]

# Plot histogram
plt.hist(lengths, bins=30)

plt.xlabel("Read length (bp)")
plt.ylabel("Number of reads")
plt.title("FASTQ Read Length Distribution")

plt.tight_layout()

# Save the figure
plt.savefig("read_length_histogram.png", dpi=300)

# Display the figure
plt.show()

------------------------------------
To create gc_content.txt
-----------------------------------
gzcat ERR10549500_1.fastq.gz | \
awk 'NR%4==2 {
    len=length($0);
    seq=$0;
    gc=gsub(/[GCgc]/,"",seq);
    print (gc/len)*100
}' > gc_content.txt


-------------------------------------------
Python code for GC content histogram
---------------------------------------
import matplotlib.pyplot as plt

with open("gc_content.txt") as f:
    gc = [float(line.strip()) for line in f]

plt.hist(gc, bins=50)

plt.xlabel("GC content (%)")
plt.ylabel("Number of reads")
plt.title("GC Content Distribution")

plt.tight_layout()
plt.savefig("gc_histogram.png", dpi=300)
plt.show()


------------------------------------------
python code for Read length with log scale 
-------------------------------
import matplotlib.pyplot as plt

# Read read lengths
with open("read_lengths.txt") as f:
    lengths = [int(line.strip()) for line in f]

# One bin for each integer read length
bins = range(min(lengths), max(lengths) + 2)

plt.figure(figsize=(10, 6))

plt.hist(
    lengths,
    bins=bins,
    align="left"
)

# Logarithmic Y-axis
plt.yscale("log")

plt.xlabel("Read length (bp)")
plt.ylabel("Number of reads (log scale)")
plt.title("FASTQ Read Length Distribution")

plt.tight_layout()

# Save figure
plt.savefig("read_length_histogram_log.png", dpi=300)

plt.show()

---------------------------------------
Python code for GC content with log scale
-------------------------------------------
import matplotlib.pyplot as plt

# Read GC percentages
with open("gc_content.txt") as f:
    gc = [float(line.strip()) for line in f]

plt.figure(figsize=(10, 6))

plt.hist(
    gc,
    bins=50
)

# Logarithmic Y-axis
plt.yscale("log")

plt.xlabel("GC content (%)")
plt.ylabel("Number of reads (log scale)")
plt.title("FASTQ GC Content Distribution")

plt.tight_layout()

# Save figure
plt.savefig("gc_histogram_log.png", dpi=300)

plt.show()









====================================================
END
====================================================