ls | grep "fq.gz$" |while read id ;do (fastqc -q -t 16 -o ./Fastqc $id );done &
java -jar /public/home/djs/software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz \
ILLUMINACLIP:/public/home/djs/software/Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:150
cd ~/huiyu/BCR_NGS/trim2
cat ../list.txt |while read id ;do (arr=($id);java -jar /public/home/djs/software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 ../${arr[0]} ../${arr[1]} -baseout ${arr[0]}.trim.fastq.gz ILLUMINACLIP:/public/home/djs/software/Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:150 );done
ls |grep "1P.fastq.gz" > 1p.txt
ls |grep "2P.fastq.gz" > 2p.txt
paste 1p.txt 2p.txt > list.txt
cat list.txt |while read id ;do arr=($id);pandaseq -f ${arr[0]} -r ${arr[1]} -F -N -o 12 -w ./${arr[0]}.merge;done
cat list.txt |while read id ;do arr=($id);pandaseq -f ${arr[0]} -r ${arr[1]} -F -N -o 12 -W ./${arr[0]}.merge;done
cat NB244-R1-H_contig.fa|grep -o CAGGCGGCC.*GAACCCAAGA |sed -n 's/CAGGCGGCC\(.*\)GAACCCAAGA/\1/p'
ls |grep "merge.out" | while read id ;do awk '{if(NR%4==2) print $0}' $id |grep -o CAGGCGGCC.*GAACCCAAGA |sed -n 's/CAGGCGGCC\(.*\)GAACCCAAGA/\1/p' > ${id}.trim_vector;done
ls |grep "trim_vector" | while read id ;do awk '{a=length($0);print $0"\t"a}' $id | sort -n -k2 |cut -f1 > ${id}.sort ;done
ls |grep "sort" | while read id ;do cat $id |uniq -c |sort -n -k1 |awk '{print ">"NR"+"$1;print $2}' > ${id}.sort.sort ;done
mkdir output_igblast
ls /public/home/djs/huiyu/BCR_NGS/trim2 |grep "sort.sort" |while read id ;do
./bin/igblastn -germline_db_V ./database/alpaca_v_igblast.fasta -germline_db_D ./database/alpaca_d_igblast.fasta -germline_db_J ./database/alpaca_j_igblast.fasta \
-query /public/home/djs/huiyu/BCR_NGS/trim2/$id -organism alpaca \
-auxiliary_data ./optional_file/alpaca_gl.aux -show_translation -extend_align5end -extend_align3end \
-num_alignments_V 1 -num_alignments_D 1 -num_alignments_J 1 \
-outfmt 19 > /public/home/djs/huiyu/BCR_NGS/trim2/output_igblast/$id.txt ;done
file=${2}
dir=${1}
/usr/bin/python3 /public/home/djs/huiyu/NGS_1.py ${dir} ${file}
cut -f 2,3,11,13,16,49 ${dir}/test1.txt > ${dir}/test2.txt
awk 'BEGIN{FS="\t"} NR > 1 {print ">"$2;print $49}' ${dir}/test1.txt > ${dir}/cluster.fa
~/software/cd-hit-v4.8.1-2019-0228/cd-hit -i ${dir}/cluster.fa -o ${dir}/cluster.test -c 0.79 -g 1 -n 4 -d 0 -l 4 -S 0 -sc 1
cat ${dir}/cluster.test.clstr |sed 's/.*>/>/g'|sed 's/\..*//g' |tr -d '\n' | sed 's/>C/\nC/g' | awk 'BEGIN{FS=">"} {for(i=2;i<=NF;i++){print $1"\t"$i}}' > ${dir}/test.cluster.txt
sed -i '1i cluster_id\tsequence_id' ${dir}/test.cluster.txt
a=$(echo ${file} |cut -d "_" -f1)
echo -e "name\tcount" >> ${file}.count.txt
echo -e "reads\t$(cat ../../Fastqc/multiqc/multiqc_data/multiqc_fastqc.txt |grep "${a}" |cut -f 5 |uniq)" >> ${file}.count.txt
ls ../ | grep ${a} |grep -v "vector" | while read id ;do echo $id && awk '{if(NR%4 == 2) print $0}' ../$id | wc -l ;done > count.txt
awk 'NR==2 {print "contig\t" $0}' count.txt >> ${file}.count.txt
echo -e "VHH_contain_contig\t$(cat test1.txt |wc -l)" >> ${file}.count.txt
cat cluster.test.clstr | grep -A 1 ">C" > represent.txt
sed -i '/--/d' represent.txt
cat represent.txt |sed 's/.*>/>/g' | sed 's/\..*//g' |awk '{if(NR%2==1) print $0}' > cluster.name.txt
cat represent.txt |sed 's/.*>/>/g' | sed 's/\..*//g' |awk '{if(NR%2==0) print $0}' > index_name.txt
paste cluster.name.txt index_name.txt > represent_cluster.txt
sed -i 's/>//g' index_name.txt
/public/home/qingw/CORE/bin/lt_create_idx cluster.fa
/public/home/qingw/CORE/bin/ltr_gwc cluster.fa index_name.txt > represent_retrive.fa
awk '{if(NR%2==1) print $0}' represent_retrive.fa > sequence_id.txt
awk '{if(NR%2==0) print $0}' represent_retrive.fa > sequence_cdr3.txt
paste sequence_id.txt sequence_cdr3.txt > represent_retrieve_1.txt
paste represent_cluster.txt represent_retrieve_1.txt | cut -f 1,4 > ${file}.dictionary.txt
echo -e "CDR3_cluster\t$(cat ${file}.dictionary.txt |wc -l)" >> ${file}.count.txt
echo -e "unique_protein\t$(cat test2.txt |cut -f5 |sort |uniq | wc -l)" >> ${file}.count.txt
/usr/bin/python3 /public/home/djs/huiyu/NGS_2.py ${dir} ${file}
import sys
import os
import numpy as np
import pandas as pd
Dir = sys.argv[1]
File = sys.argv[2]
os.chdir(Dir)
df1 = pd.read_csv(File,sep='\t')
df2 = df1[(df1["stop_codon"] == 'F') & (df1["vj_in_frame"] == 'T') & (df1["v_frameshift"] == 'F') & (df1["productive"] == 'T') & (df1["complete_vdj"] == 'T')]
df3 = df2.sort_values(["v_call","j_call"])
df3.to_csv("test1.txt",sep="\t")
import sys
import os
import numpy as np
import pandas as pd
Dir = sys.argv[1]
File = sys.argv[2]
os.chdir(Dir)
df1 = pd.read_csv("test2.txt",sep="\t")
df2 = pd.read_csv("test.cluster.txt",sep="\t")
df3 = pd.merge(df1,df2,on="sequence_id")
df3["cDNA_contig"] = df3["sequence_id"].map(lambda x:x.split('+')[1])
df3.cDNA_contig = df3.cDNA_contig.astype("int64")
df4 = df3.groupby(["v_call","j_call","cluster_id"]).cDNA_contig.sum()
df4.name = "counts_of_contig_in_the_same_clonotype"
df5 = pd.merge(df3,df4,on=["v_call","j_call","cluster_id"])
df5["clonotype"] = df5["cluster_id"].map(lambda x:x.split(" ")[1])+"@"+df5["v_call"].map(lambda x:x.split("*")[0])+"@"+df5["j_call"].map(lambda x:x.split("*")[0])
df6 = df5.groupby("cluster_id").cDNA_contig.sum()
df6.name = "counts_of_contig_in_the_same_cluster"
df6 = pd.merge(df5,df6,on="cluster_id")
df7= df3.groupby("sequence_alignment_aa").cDNA_contig.sum()
df7.name = "counts_of_contig_in_the_same_AA"
df7 = pd.merge(df7,df6,on="sequence_alignment_aa")
df7 = df7.loc[:,["clonotype","counts_of_contig_in_the_same_clonotype","cluster_id","counts_of_contig_in_the_same_cluster","cdr3_aa","counts_of_contig_in_the_same_AA","cDNA_contig","v_call","j_call","sequence_alignment_aa","sequence"]]
df7.cluster_id = df7.cluster_id.map(lambda x:x.split(' ')[1])
df7["cluster_id"] = df7.cluster_id.astype("int64")
df7["cluster_id"] = df7.cluster_id + 1
df8 = df7.sort_values(["counts_of_contig_in_the_same_AA","cDNA_contig"],ascending=False)
df8 = df8.drop_duplicates(subset="sequence_alignment_aa",keep="first")
df9 = df8.loc[df8.counts_of_contig_in_the_same_AA >= 10]
df10 = df9.sort_values(["cDNA_contig","counts_of_contig_in_the_same_AA","counts_of_contig_in_the_same_cluster"],ascending=False)
df10 = df10.drop_duplicates(subset="clonotype", keep='first')
df11 = df9.sort_values(["counts_of_contig_in_the_same_cluster","cDNA_contig","counts_of_contig_in_the_same_AA"],ascending=False)
df11 = df11.drop_duplicates(subset="cluster_id", keep='first')
df12 = df11.loc[:,["cluster_id","counts_of_contig_in_the_same_cluster","cdr3_aa"]]
df12["frequency_of_total"] = df12.counts_of_contig_in_the_same_cluster/df3.cDNA_contig.sum()
df12["frequency"] = df12.counts_of_contig_in_the_same_cluster/df12.counts_of_contig_in_the_same_cluster.sum()
df12["CDR3_length"] = df12.cdr3_aa.apply(lambda x:len(x))
df12 = df12.loc[:,["cluster_id","counts_of_contig_in_the_same_cluster","frequency","frequency_of_total","cdr3_aa","CDR3_length"]]
df12 = df12.sort_values("cluster_id")
df13 = pd.read_csv(File+".dictionary.txt",sep="\t",names=["cluster_id","represent"])
df13["cluster_id"] = df13["cluster_id"].map(lambda x:x.split(' ')[1])
df13.cluster_id = df13.cluster_id.astype("int64")
df13.cluster_id = df13.cluster_id + 1
df12.cluster_id = df12.cluster_id.astype("string")
df13.cluster_id = df13.cluster_id.astype("string")
df14 = pd.merge(df12,df13,on="cluster_id",how="inner")
df16 = df8.groupby("cluster_id",as_index=False).agg("count")
df16 = df16.loc[:,["cluster_id","sequence_alignment_aa"]]
df16.rename(columns={"sequence_alignment_aa":"variety_of_protein"})
df16.cluster_id = df16.cluster_id.astype("string")
df14 = pd.merge(df14,df16,on="cluster_id",how="inner")
df15 = pd.read_csv(File+".count.txt",sep="\t")
df15.iloc[2,1] = df3.cDNA_contig.sum()
with pd.ExcelWriter(File+".result.xlsx") as writer:
df7.to_excel(writer, sheet_name="all_of_contig",index=False)
df8.to_excel(writer, sheet_name="table_of_unique_protein",index=False)
df10.to_excel(writer, sheet_name="table_of_unique_clonotype",index=False)
df11.to_excel(writer, sheet_name="table_of_unique_cluster",index=False)
df14.to_excel(writer, sheet_name="table_of_analysis",index=False,startrow=8)
df15.to_excel(writer,sheet_name="table_of_analysis",index=False,startrow=0)