> For the complete documentation index, see [llms.txt](https://zpliu.gitbook.io/booknote/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/shu-ju-chu-li-liu-cheng/20200323.md).

# 第三个结果

## 统计各个棉种中剪切数据

```bash
awk '$3~/IntronR/{a+=1}$3~/ExonS/{b+=1}$3~/AltA/{c+=1}$3~/AltD/{d+=1}$3~/AltP/{e+=1}$3~/Other/{f+=1}END{print a,b,c,d,e,f}' OFS="\t"  end_third
```

> TM-1去除了scaffold

| 棉种   | IR    | ExonS | AltA  | AltD  | AltP | Other | Total |
| ---- | ----- | ----- | ----- | ----- | ---- | ----- | ----- |
| A2   | 28078 | 2276  | 5199  | 4135  | 1823 | 1435  | 42946 |
| D5   | 31756 | 4280  | 7055  | 4900  | 1819 | 1805  | 51615 |
| At   | 26763 | 2328  | 6063  | 5188  | 2599 | 1448  |       |
| Dt   | 27412 | 2431  | 6227  | 5371  | 2658 | 1625  |       |
| TM-1 | 54175 | 4759  | 12290 | 10559 | 5257 | 3073  | 90113 |

对应的基因数

```bash
for i in IntronR ExonS AltA AltD AltP Other
do
awk '$3=="'"$i"'"{print $2}' end_third|sort |uniq|wc -l
done
```

| 棉种   | IR    | ExonS | AltA | AltD | AltP | Other | total |
| ---- | ----- | ----- | ---- | ---- | ---- | ----- | ----- |
| A2   | 9203  | 1560  | 3376 | 2810 | 1296 | 949   |       |
| D5   | 10581 | 2958  | 4731 | 3427 | 1303 | 1209  |       |
| At   | 8901  | 1679  | 4000 | 3480 | 1701 | 989   |       |
| Dt   | 9143  | 1723  | 4092 | 3557 | 1712 | 1044  |       |
| TM-1 | 18044 | 3402  | 8092 | 7037 | 3413 | 2033  |       |

## 统计每个isform发生IR的次数

## 提取各个棉种中同源基因的剪切事件

```bash
cut -f4 ../A2_vs_D5_collinearity.txt|xargs -I {} grep {} ../../D5/end_third|awk '{print $1,$4,$5,$2"-"$1"-"$4"-"$5,$3}' OFS="\t"
```

从文件夹中*static.txt*文件中提取

```bash
```

### 提取独有的基因对应的剪切事件

```bash
awk '{print $1"\t"$3}' At_vs_Dt_collinerity.txt|cat - ~/work/Alternative/result/duplicate_gene/At_Dt/At_vs_Dt.gff|awk '$2~/Ghir_A/{print $2}' |sort |uniq -u
```

### 进行Blast保守性分析

提取左右两端300bp序列

```bash
#提取5'端序列
awk -F "\t" '$4=="+"{print $1,$2-300,$2-1,$4,$5}$4=="-"{print $1,$3+1,$3+300,$4,$5}' OFS="\t" ../D5_IR.txt >5_sequence.bed
#提取3’端序列
awk -F "\t" '$4=="-"{print $1,$2-300,$2-1,$4,$5}$4=="+"{print $1,$3+1,$3+300,$4,$5}' OFS="\t" ../D5_IR.txt >3_sequence.bed
```

合并两端序列 脚本位置 `~/script/Alternative/merge.py`

进行`All-vs-All blast`

```bash
file="A2_D5"
echo "blastn -query blastDB/${file}.fasta -db blastDB/${file} -evalue 1e-5 -num_threads 10 -outfmt 6 -out ${file}.blast"
file="A2_D5"
bsub -J Dt_Dt -q "normal" -n 10 -R span[hosts=1] -e Dt_Dt.err -o Dt_Dt.out "`echo blastn -query blastDB/${file}.fasta -db blastDB/${file} -evalue 1e-5 -num_threads 10 -outfmt \'6  qseqid sseqid qstart qend sstart send nident pident qcovs evalue bitscore\' -out ${file}.blast`"
```

### 提取保守的IR

* 匹配长度大于200bp
* 相似度大于90%
* 覆盖度超过85%

```bash
awk '$7>200&&$8>90&&$9>85{print $0}'
```

### 使用intron比对到gene上

1.提取对应的同源基因序列

```bash
cut -f3 A2_vs_At_collinearity.txt|xargs -I {} grep {} ~/work/Alternative/data/Ga_genome/G.arboreum.Chr
.v1.0.gff|awk '$3~/gene/{print $1,$4,$5,$7,$9}' OFS="\t" |awk -F ";" '{print $1}' |sed 's/ID=//g' |awk '{print $1,$2,$3,$4,"gene-"$5"-"$1"-"$2"-"$3"-"$4}' OFS="\t" >A2_gene.bed
```

2.提取对应的IR序列

3.将intron序列比对到基因区域，进行筛选

* 相似度90以上
* 覆盖度90以上
* 首先是同源基因

```bash
python /public/home/zpliu/work/Alternative/result/duplicate_gene/IR_align_to_gene.py -homologous ../A2_vs_At_collinearity.txt -Blast 1 -out 2
#提取比对到基因区域的IR
awk '$4~/evm/{print $1,$2,$3,"-",$4}' OFS="\t" 2 >A2_IR_to_Atgene.bed
awk '$4~/Ghir/{print $1,$2,$3,"-",$4}' OFS="\t" 2 >At_IR_to_A2gene.bed
```

**这里发现有些IR的intron不能够比对到对应的同源基因上，不知道是不是筛选指标太严格了**

* 把指标松一些，比对上的片段占intron的0.8即可。
* 分为比对上的intron
  * 对应同源基因有IR
  * 对应同源基因没有IR

```bash
#与对应同源基因的IR取交集
~/software/bedtools2-2.29.0/bin/intersectBed -loj -a A2_IR_to_Atgene.bed -b ../At_IR.txt |cut -f5,10 >3
#类似
~/software/bedtools2-2.29.0/bin/intersectBed -loj -a At_IR_to_A2gene.bed -b ../A2_IR.txt|cut -f5,10 >4
#合并两个文件
awk '$2!="."{print $2"\t"$1}$2=="."{print $0}' 4|cat - 3|sort |uniq >A2_At_conserveIR.txt 

evm.TU.Ga01G0020-Chr01-129583-130080--    .
evm.TU.Ga01G0021-Chr01-142321-143057-+    .
evm.TU.Ga01G0021-Chr01-142321-143808-+    Ghir_A01G000210-Ghir_A01-157898-158012-+
evm.TU.Ga01G0021-Chr01-142321-143808-+    Ghir_A01G000210-Ghir_A01-158029-158264-+
Ghir_A13G024820-Ghir_A13-108246222-108246433--    .
Ghir_A13G024830-Ghir_A13-108259086-108259204-+    .
```

> 其中点表示，对应同源基因序列类似，但是附近没有IR事件
>
> 不是点则代表保守的intronR

* 没有比对上的intron

  > 感觉应该也不会保守，就不管了；直接统计出来

```bash
cut -f5 A2_IR_to_Atgene.bed|sort |uniq >666
cut -f5 ../A2_IR.txt|sort |uniq  >777
cat 666 777|sort |uniq -u  >A2_unmapping_IR
cut -f5 At_IR_to_A2_gene.bed|sort |uniq >666
cut -f5 ../At_IR.txt|sort |uniq >777
cat 666 777|sort |uniq -u >At_unmapping_IR
```

### 整套流程

```bash
file1=At
file2=Dt
cat ../${file1}_gene.bed.fa ../${file2}_gene.bed.fa >${file1}_${file2}_gene.fa
cat ../${file1}_IR.txt.fa  ../${file2}_IR.txt.fa >${file1}_${file2}_IR.fa
makeblastdb -dbtype "nucl" -in ${file1}_${file2}_gene.fa -out ${file1}_${file2}
blastn -query ${file1}_${file2}_IR.fa -db ${file1}_${file2} -evalue 1e-5 -num_threads 10 -outfmt "6  qseqid sseqid qstart qend sstart send nident pident qcovs evalue bitscore" -out 1
python /public/home/zpliu/work/Alternative/result/duplicate_gene/IR_align_to_gene.py -homologous ../${file1}_vs_${file2}_collinearity.txt -Blast 1 -out 2

##这里要个性化一些，哪个命名在前，awk中就用哪个基因编号
awk '$4~/evm/{print $1,$2,$3,"-",$4}' OFS="\t" 2 >${file1}_IR_to_${file2}gene.bed
awk '$4~/Ghir/{print $1,$2,$3,"-",$4}' OFS="\t" 2 >${file2}_IR_to_${file1}gene.bed
###
~/software/bedtools2-2.29.0/bin/intersectBed -loj -a ${file1}_IR_to_${file2}gene.bed -b ../${file2}_IR.txt |cut -f5,10 >3
~/software/bedtools2-2.29.0/bin/intersectBed -loj -a ${file2}_IR_to_${file1}gene.bed -b ../${file1}_IR.txt|cut -f5,10 >4
awk '$2!="."{print $2"\t"$1}$2=="."{print $0}' 4|cat - 3|sort |uniq >${file1}_${file2}_conserveIR.txt 
cut -f5  ${file1}_IR_to_${file2}gene.bed|sort |uniq >666
cut -f5 ../${file1}_IR.txt|sort |uniq  >777
cat 666 777|sort |uniq -u  >${file1}_unmapping_IR
cut -f5  ${file2}_IR_to_${file1}gene.bed|sort |uniq >666
cut -f5 ../${file2}_IR.txt|sort |uniq >777
cat 666 777|sort |uniq -u >${file2}_unmapping_IR
```