In my continuing quest to conquer fastQ, fastA, sam & bam files, I have accumulated several useful ‘tools’. Many of them are included in other software packages (e.g. SAMtools), but for some tasks, especially file management and  conversion, no standard toolkit exists, and instead researchers script their own solution.

For me, sed and awk, along with a few other standard *nix tools have been extremely useful, a few of them are generally useful, I think… In hopes of helping others, I’m going to start a list, here. The initial list is short (these were the ones I used in the past week or so, fresh in my mind), but I’ll plan to add and update as I discover new, better ways of handling these files. MOst all of these were inspired by stuff I’ve seen on the internet, so if you have something you want to share.. please add it in the comments!

#Create an unwrapped fasta file, you should add ambiguous nucleotides to the search if you have them, also, fasta defline ending in ‘ACTGNn-‘ will trip the script up, so use caution…

sed -i ':begin;$!N;/[ACTGNn-]n[ACTGNn-]/s/n//;tbegin;P;D' test.fa

#Add /1 to end of fastQ defline, good for assemblers that need those to tell left and right files. Note that you should change [0-9] to whatever follows your @– this is machine specific.

sed -i 's_^@[0-9]:.*_&/1_g' left.fq
sed -i 's_^@[0-9]:.*_&/2_g' right.fq

#Deinterleave a shuffled fastq file

sed -n '/2:N/{N;N;N;p;}' shuff.fastq > right.fq
sed -n '/1:N/{N;N;N;p;}' shuff.fastq > left.fq

##Get the sum of a column (Column 2 ($2)) in this case

awk '{add +=$2} END {print add}'

#Get SD of a column (Column 2 ($2) in this case

awk 'NR>2 {sum+=$2; array[NR]=$2} END {for(x=1;x<=NR;x++){sumsq+=((array[x]-(sum/NR))^2);}print sqrt(sumsq/NR)}'

#Print mean of column (Column 2 ($2)) in this case

awk '{sum+=$2} END { print "Average = ",sum/NR}'

#Remove duplicate entries in column 10, and print new spreadsheet at unique10.txt

cat test.txt | sort -k10 | awk '!a[$10]++' > unique10.txt

#Retain entries in speadsheet where column 3 is greater than 300 (imagine, culling by sequence length)

cat table.txt | awk '300>$3{next}1' > morethan300.txt

#Extract fastQ from BAM for a bunch of BAM files– requires Picard. (sorry, no sed or awk in this one)

for i in `ls *bam`;
do F=`basename $i .bam`;
java -Xmx9g -jar ~/software/SamFormatConverter.jar MAX_RECORDS_IN_RAM=1000000 INPUT=$i OUTPUT=/dev/stdout | java -Xmx9g -jar ~/software/SamToFastq.jar MAX_RECORDS_IN_RAM=1000000 INPUT=/dev/stdin FASTQ=$F.1.fq SECOND_END_FASTQ=$F.2.fq;
done