## 1 解压缩
tar -xzf 03.tgz

## 2. 从03/文件夹里获取指定物种的序列

for i in 03/*_nucl.fasta
  do
        exon_id=`basename ${i%_nucl.fasta}`
        seq=`cat ${i} | grep -A 1 "metzeb" | tail -n 1 | sed 's/\-//g'`
        seq_length=`echo ${seq} | wc -m` ###计算序列的非-的长度，
        # Include sequence records only if the sequence is not completely missing.
        if [[ ${seq_length} -gt 1 ]]
        then
                echo ">${exon_id}"
                echo "${seq}"
        fi
done > metzeb.fasta 

## 3. 对指定物种的fastq.gz解压缩 
gunzip neopul.r1.fastq.gz
gunzip neopul.r2.fastq.gz

## 4. astram第一步程序 generate read library in SQL format 
### -b prefix
python3 aTRAM/atram_preprocessor.py -b neopul --end-1 neopul.r1.fastq --end-2 neopul.r2.fastq

## 5. astram第二步程序 对指定的metzeb.fasta序列进行定向组装
### -a 选择abyss/ trinity/ spides/ velvet的一个
python3 aTRAM/atram.py -b neopul -Q metzeb.fasta -a abyss -o neopul

## 6. 将找到的所有orthologs输出到一个文件
cat neopul.neopul_metzeb_*.filtered_contigs.fasta > neopul.fasta

## 1. 建库
makeblastdb -in neopul.fasta -dbtype nucl

## 2. 利用zebra fish的序列检测
python3 find_orthologs.py -s 1 metzeb.fasta neopul.fasta

##3. 后续
ls ENSDARE*.fasta | wc -l

for i in ENSDARE*.fasta
  do
        tail -n 1 ${i} | grep -e A -e C -e G -e T
done | wc -l

 mkdir 04
  for i in 03/*_nucl.fasta
  do
  	exon_id=`basename ${i%_nucl.fasta}`
  	if [ -f ${exon_id}.fasta ]
  	then
  		# Copy only the neopul sequence to a temporary file.
  		tail -n 2 ${exon_id}.fasta > tmp.fasta
  		# Add the neopul sequence to the existing alignment of 11 cichlids.
  		mafft --inputorder --keeplength --anysymbol --add tmp.fasta ${i} > 04/${exon_id}_nucl.fasta 
  		# Clean up the temporary file.
  		rm tmp.fasta
  	fi
done


## 2. 去除neo类群的序列，为了后面的Bayesian analysis of species networks
mkdir 05
  for i in 04/*_nucl.fasta
  do
        exon_id=`basename ${i%_nucl.fasta}`
        python3 convert.py ${i} 05/${exon_id}_nucl.fasta -f fasta -p neo
done

## 3. 利用filter_exons_by_missing_data.rb 过滤
ruby filter_exons_by_missing_data.rb 05 06 0 0

## 4. 利用BMGE去除低质量比对区域
ruby filter_sites_with_BMGE.rb 06 07 BMGE.jar 0.2 0.5

## 5. 去除gap区域
ruby filter_exons_by_missing_data.rb 07 08 0 150

## 6. 合并同属一个gene的exon序列
ruby concatenate_exons_per_gene.rb 08 09 exons_info.txt

## 7. 保留500 bp以上序列

mkdir 10
  for i in 09/*.fasta
  do
        gene_id=`basename ${i%.fasta}`
        seq_length=`head -n 2 ${i} | tail -n 1 | wc -m`
        if [[ ${seq_length} -ge 500 ]]
        then
                cp ${i} 10/${gene_id}.fasta
        fi
done