植物长链非编码RNA（lncRNA）鉴定实例（拟南芥数据）

SEQLIBS=(SRR8428909 SRR8428908 SRR8428906 SRR8428905 )
for seqlib in ${SEQLIBS[@]}; do    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR842/00${seqlib:9:10}/${seqlib}/${seqlib}_1.fastq.gz    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR842/00${seqlib:9:10}/${seqlib}/${seqlib}_2.fastq.gzdone

bash download_raw_data.sh

bgzip -d SRR8428909_1.fastq.gzbgzip -d SRR8428909_2.fastq.gz
mv SRR8428909_1.fastq wt_Rep1_R1.fastqmv SRR8428909_2.fastq wt_Rep1_R2.fastq
bgzip -d SRR8428908_1.fastq.gzbgzip -d SRR8428908_2.fastq.gz
mv SRR8428908_1.fastq wt_Rep2_R1.fastqmv SRR8428908_2.fastq wt_Rep2_R2.fastq
bgzip -d SRR8428906_1.fastq.gzbgzip -d SRR8428906_2.fastq.gz
mv SRR8428906_1.fastq EE_Rep1_R1.fastqmv SRR8428906_2.fastq EE_Rep1_R2.fastq
bgzip -d SRR8428905_1.fastq.gzbgzip -d SRR8428905_2.fastq.gz
mv SRR8428905_1.fastq EE_Rep2_R1.fastqmv SRR8428905_2.fastq EE_Rep2_R2.fastq

SEQLIBS=(EE_Rep1 EE_Rep2  wt_Rep1 wt_Rep2 )
for seqlib in ${SEQLIBS[@]}; do    fastp -i ${seqlib}_R1.fastq -I ${seqlib}_R2.fastq -o ${seqlib}_clean_R1.fastq -O ${seqlib}_clean_R2.fastqdone

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gzbgzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gzmv  Arabidopsis_thaliana.TAIR10.dna.toplevel.fa At.fawget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.40.gff3.gzbgzip -d Arabidopsis_thaliana.TAIR10.40.gff3.gz mv Arabidopsis_thaliana.TAIR10.40.gff3 At.gff3

mkdir referencemv At* referencecd referencebowtie2-build At.fa At

cd ../tophat2 -p 8 -I 5000 -G reference/At.gff3 -o wt1_thout reference/At ../wt_Rep1_clean_R1.fastq ../wt_Rep1_clean_R2.fastqtophat2 -p 12 -I 5000 -G reference/At.gff3 -o wt2_thout reference/At ../wt_Rep2_clean_R1.fastq ../wt_Rep2_clean_R2.fastq
tophat2 -p 8 -I 5000 -G reference/At.gff3 -o EE1_thout reference/At ../EE_Rep1_clean_R1.fastq ../EE_Rep1_clean_R2.fastqtophat2 -p 8 -I 5000 -G reference/At.gff3 -o EE2_thout reference/At ../EE_Rep2_clean_R1.fastq ../EE_Rep2_clean_R2.fastq

cufflinks -p 4 -g reference/At.gff3 -I 5000 -o wt1_clout wt1_thout/accepted_hits.bamcufflinks -p 4 -g reference/At.gff3 -I 5000 -o wt2_clout wt2_thout/accepted_hits.bamcufflinks -p 4 -g reference/At.gff3 -I 5000 -o EE1_clout EE1_thout/accepted_hits.bamcufflinks -p 4 -g reference/At.gff3 -I 5000 -o EE2_clout EE2_thout/accepted_hits.bam

find . -name transcripts.gtf > assemblies.txtcuffmerge -p 4 -g reference/At.gff3 -s reference/At.fa assemblies.txt

cuffdiff -o diff_out -b reference/At.fa -p 8 -T -L EE,wt -u merged_asm/merged.gtf EE1_thout/accepted_hits.bam,EE2_thout/accepted_hits.bam wt1_thout/accepted_hits.bam,wt2_thout/accepted_hits.bam

cat merged_asm/merged.gtf | grep 'class_code "[uiox]"' > selected.gtfgffread -w selected.fa -g reference/At.fa selected.gtfcat selected.fa | grep ">" | wc -l

import sysfrom Bio import SeqIO
inputfasta = sys.argv[1]min_len = int(sys.argv[2])outputfasta = sys.argv[3]
i = 0j = 0
fw = open(outputfasta,'w')
for rec in SeqIO.parse(inputfasta,'fasta'):	i += 1	if len(rec.seq) > min_len:		j += 1		fw.write(">%s\n%s\n"%(rec.id,str(rec.seq)))		print("The total number of sequence is: ",i)print("Number of sequence retained is: ",j)

python .\select_seq_according_to_seq_length.py .\selected.fa 200 output.fasta

cat result_cpc2.txt | grep 'noncoding' | cut -f 1 > non-coding-transcript.txtwc -l non-coding-transcript.txt