Commit 61d48924 authored by jplantad's avatar jplantad
Browse files

get_consensus_sequence_bygene.sh and align_with_mafft_bygene.sh

parent 7ec48582
#!/bin/bash
# align human CDS with BST2 sequences
# set the humain CDS path and file
path_hgCDS="/home/stagiaire/Bureau/data_lucie/apeVIP_Lucie_02052021/to_get_cds_fasta_seq/genomic_DNA_sequences/CDS/"
hg_CDS="BST2_hg18_CDS.fasta"
# set the sequences path and file
path_seq="/home/stagiaire/Bureau/phylogenetics/data_sequences/pon/"
seq="BST2_ponAbe2_all_RC.fas"
# set the output and merge sequences with human CDS
path_out="/home/stagiaire/Bureau/phylogenetics/data_sequences/mafft_alignment/"
out_cat="BST2_ponAbe2-hg18_all_RC_CDS.fa"
cat ${path_seq}${seq} ${path_hgCDS}${hg_CDS} > ${path_out}${out_cat}
# align sequences with human CDS using mafft
out_mafft="BST2_ponAbe2-hg18_RC_CDSmafft.fa"
mafft ${path_out}${out_cat} > ${path_out}${out_mafft}
#!/bin/bash
# get 2 consensus sequence per individual using bcftools
# help : https://samtools.github.io/bcftools/bcftools.html#consensus
home="/home/stagiaire/Bureau/phylogenetics/"
# set the vcf.gz file
path_out=${home}"data_sequences/pon/"
cd ${home}
file_vcf="999_BST2_ponAbe2.vcf_fixed.vcf.gz"
# set coordinates
coord="chr19:17793264-17795651"
# list individuals
list_indiv=$(cat ${home}"indiv_list/pon_indiv_list.txt")
for indiv in ${list_indiv}
do
tmp=$(echo ${indiv} | tr "/" "\n")
tmp2=$(echo "${tmp%%.*}" | sed -n 2p)
header="BST2_"${tmp2}
# get the first sequence
samtools faidx ${home}ref_genomes/ochPri2.agp ${coord} | bcftools consensus -H1 -s ${indiv} ${path_out}${file_vcf} -o ${path_out}${header}"_seqA_alt.fa"
# replace the header with individual info
sed "s/>/>${tmp2}_A_/g" ${path_out}${header}"_seqA_alt.fa" > ${path_out}${header}"_renamed_AB_alt.fa"
# get the second sequence
samtools faidx ${home}ref_genomes/ochPri2.agp ${coord} | bcftools consensus -H2 -s ${indiv} ${path_out}${file_vcf} -o ${path_out}${header}"_seqB_alt.fa"
# replace the header with individual info, and merge it with seqA
sed "s/>/>${tmp2}_B_/g" ${path_out}${header}"_seqB_alt.fa" >> ${path_out}${header}"_renamed_AB_alt.fa"
done
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment